{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. 导入工具包 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import scipy.io as sio\n",
    "import scipy.sparse as ss\n",
    "\n",
    "from sklearn.preprocessing import normalize\n",
    "from sklearn.cluster import MiniBatchKMeans\n",
    "from sklearn import metrics\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. 读入数据"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 读入PE_eventIndex.pkl，获取train.csv和test.csv中出现的events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpath = \"./data/\"\n",
    "#train = pd.read_csv(dpath + \"train.csv\")\n",
    "#test = pd.read_csv(dpath + \"test.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train events count: 8846\n",
      "Test events count: 6173\n",
      "Total events count in train and test: 15019\n",
      "Unique events count in train & test: 13418\n"
     ]
    }
   ],
   "source": [
    "# 用set获取unique events\n",
    "#uniqueEvents = set()\n",
    "\n",
    "#train_values = set(train[\"event\"].values)\n",
    "#test_values = set(test[\"event\"].values)\n",
    "\n",
    "#print(\"Train events count: %d\" % len(train_values))\n",
    "#print(\"Test events count: %d\" % len(test_values))\n",
    "#print(\"Total events count in train and test: %d\" % (len(train_values) + len(test_values)))\n",
    "\n",
    "#uniqueEvents = train_values | test_values\n",
    "#print(\"Unique events count in train & test: %d\" % len(uniqueEvents))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total events count in train & test :13418\n"
     ]
    }
   ],
   "source": [
    "# 读取Users_Events生成的字典内容\n",
    "eventIndex = pickle.load(open(\"PE_eventIndex.pkl\", 'rb'))\n",
    "n_events = len(eventIndex)\n",
    "\n",
    "print(\"Total events count in train & test :%d\" % n_events)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 逐行读入events.csv，将train和test中出现过的event相关内容取出"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 打开文件\n",
    "events_file = open(dpath + \"events.csv\", \"r\")\n",
    "\n",
    "# 读取第一行的列名\n",
    "col_names = events_file.readline().strip()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 生成稀疏矩阵\n",
    "eventContMatrix = ss.dok_matrix((n_events, 101))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 逐行读入events.csv\n",
    "for line in events_file.readlines():\n",
    "    cols = line.strip().split(\",\")\n",
    "    eventId = str(cols[0])\n",
    "    \n",
    "    # 判断是否出现在train或test中\n",
    "    if eventId in eventIndex:\n",
    "        i = eventIndex[eventId]\n",
    "        \n",
    "        for j in range(9, 110):\n",
    "            eventContMatrix[i, j-9] = cols[j]\n",
    "\n",
    "events_file.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 3.,  1.,  2., ...,  0.,  0., 54.],\n",
       "       [ 4.,  3.,  2., ...,  0.,  0., 37.],\n",
       "       [ 0.,  1.,  0., ...,  0.,  0.,  6.],\n",
       "       ...,\n",
       "       [ 1.,  0.,  0., ...,  0.,  0., 23.],\n",
       "       [ 1.,  0.,  0., ...,  0.,  0., 36.],\n",
       "       [ 2.,  2.,  0., ...,  0.,  0., 76.]])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eventContMatrix.toarray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 用L2模归一化\n",
    "eventContMatrix = normalize(eventContMatrix, norm = \"l2\", axis = 0, copy = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.00132991, 0.00261414, 0.00570798, ..., 0.        , 0.        ,\n",
       "        0.00373071],\n",
       "       [0.00177322, 0.00784242, 0.00570798, ..., 0.        , 0.        ,\n",
       "        0.00255622],\n",
       "       [0.        , 0.00261414, 0.        , ..., 0.        , 0.        ,\n",
       "        0.00041452],\n",
       "       ...,\n",
       "       [0.0004433 , 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.001589  ],\n",
       "       [0.0004433 , 0.        , 0.        , ..., 0.        , 0.        ,\n",
       "        0.00248714],\n",
       "       [0.00088661, 0.00522828, 0.        , ..., 0.        , 0.        ,\n",
       "        0.00525062]])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eventContMatrix.toarray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "sio.mmwrite(\"EV_eventContMatrix\", eventContMatrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. 对活动进行聚类"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 进行类别数目为K的聚类模型，并返回评价结果\n",
    "def K_cluster_analysis(K, df):\n",
    "    print(\"K-means begin with %d clusters\" % K)\n",
    "    \n",
    "    # K-means 训练\n",
    "    km = MiniBatchKMeans(n_clusters = K)\n",
    "    km.fit(df)\n",
    "    \n",
    "    # 保存预测结果\n",
    "    cluster_result = km.predict(df)\n",
    "    \n",
    "    # 评估结果\n",
    "    #CH_score = metrics.calinski_harabaz_score(df, cluster_result)\n",
    "    CH_score = metrics.silhouette_score(df,cluster_result)\n",
    "    print(\"CH_score: {}\".format(CH_score))\n",
    "    \n",
    "    return CH_score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "calinski_harabaz_score计算速度比silhouette_score快很多；calinski_harabaz取得的值很大，而silhouette_score在-1 到 1 之间；两种评价指标的大致趋势一致。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K-means begin with 10 clusters\n",
      "CH_score: 0.013842842471957197\n",
      "K-means begin with 20 clusters\n",
      "CH_score: -0.1035501559489174\n",
      "K-means begin with 30 clusters\n",
      "CH_score: -0.12148696349479626\n",
      "K-means begin with 40 clusters\n",
      "CH_score: -0.04654268191903653\n",
      "K-means begin with 50 clusters\n",
      "CH_score: -0.07304087418477744\n",
      "K-means begin with 60 clusters\n",
      "CH_score: -0.0919474539887508\n",
      "K-means begin with 70 clusters\n",
      "CH_score: -0.0498539891522519\n",
      "K-means begin with 80 clusters\n",
      "CH_score: -0.04348239684620829\n",
      "K-means begin with 90 clusters\n",
      "CH_score: -0.18931584665717002\n",
      "K-means begin with 100 clusters\n",
      "CH_score: -0.15040276864393018\n"
     ]
    }
   ],
   "source": [
    "# 设置超参数（聚类数目K）的搜索范围\n",
    "CH_scores = []\n",
    "Ks = [10,20,30,40,50,60,70,80,90,100]\n",
    "for K in Ks:\n",
    "    ch = K_cluster_analysis(K, eventContMatrix.toarray())\n",
    "    CH_scores.append(ch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.013842842471957197, -0.1035501559489174, -0.12148696349479626, -0.04654268191903653, -0.07304087418477744, -0.0919474539887508, -0.0498539891522519, -0.04348239684620829, -0.18931584665717002, -0.15040276864393018]\n"
     ]
    }
   ],
   "source": [
    "print(CH_scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x3300944b00>]"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD8CAYAAABkbJM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XuclHX5//HXhQsKYgKCiIqIhmfDw4qijpkIikyiaab5VSzJjl/Ks2Vl38xfipWmZqVoYpknPICEB8RKPLuogIQKiiaKuMpBPCJ4/f647o2VZtnZncM9u/N+Ph7zmLnvuee+r52dnWs/Z3N3REREmtMh7QBERKRtUMIQEZG8KGGIiEhelDBERCQvShgiIpIXJQwREcmLEoaIiORFCUNERPKihCEiInmpSTuAYurZs6dvvfXWaYchItKmzJgx4y1379Xcce0qYWy99dbU1dWlHYaISJtiZq/kc5yqpEREJC9KGCIikhclDBERyYsShoiI5EUJQ0RE8qKEISIieVHCEBGRvChhAHPmwGmnwYcfph2JiEjlUsIAXn4ZLrkE/vnPtCMREalcShjAQQdB587wt7+lHYmISOVSwiCSxZAhMHkyuKcdjYhIZVLCSGSzsGABzJ2bdiQiIpVJCSNx2GFxP3lyunGIiFQqJYxE374wcKAShohIU5QwGslm4ZFHYMmStCMREak8ShiNZLOwejXce2/akYiIVB4ljEb22gt69lS1lIhILkoYjay3XjR+3303rFqVdjQiIpWlaAnDzA41s+fNbL6ZnZPj+fXN7Obk+cfNbOtGz/0w2f+8mR2S7zlLIZuFpUvh0UfLcTURkbajKAnDzNYDfgcMB3YCjjOzndY67GRgqbt/FrgEuCh57U7AscDOwKHAlWa2Xp7nLLphw6CmRqO+RUTWVqwSxiBgvru/5O4rgZuAkWsdMxIYnzyeAAwxM0v23+TuH7n7AmB+cr58zll0G28MBxygdgwRkbUVK2FsAbzaaHthsi/nMe6+ClgObLKO1+ZzzpIYMSJmsF2woBxXExFpG4qVMCzHvrVnZWrqmJbu//RJzU4xszozq6uvr2820Hxks3GvaikRkTWKlTAWAn0bbW8JvN7UMWZWA2wMLFnHa/M5J+5+lbvXunttr169CvwxwnbbwYABShgiIo0VK2E8CQwws/5m1oloxJ601jGTgFHJ46OBB9zdk/3HJr2o+gMDgCfyPGfJZLPwwAPw7rvluqKISGUrSsJI2iS+B9wLzAVucfc5ZvZzMzs8OewaYBMzmw+cBpyTvHYOcAvwL+Ae4LvuvrqpcxYj3nxks7ByJUybVq4riohUNvN2tABEbW2t19XVFeVcK1fGqO+vfAWuvroopxQRqUhmNsPda5s7TiO9m9CpExxySLRjtKOcKiLSakoY65DNwqJF8PTTaUciIpI+JYx1GD4czDSIT0QElDDWadNNYe+9lTBEREAJo1kjRsCTT8Ibb6QdiYhIupQwmtEw6nvKlHTjEBFJmxJGMwYOhC220KhvEREljGaYRSnjvvvgo4/SjkZEJD1KGHnIZmOKkAcfTDsSEZH0KGHk4aCDYIMN1FtKRKqbEkYeunSJpDF5skZ9i0j1UsLIUzYLL70Ezz+fdiQiIulQwsjTiBFxr2opEalWShh52mor+NznlDBEpHopYbTAiBHw0EOwdGnakYiIlF9BCcPMepjZVDObl9x3b+K4Uckx88xsVLKvi5n9zcyeM7M5ZnZho+NPMrN6M3smuY0uJM5iyWZh9eoYkyEiUm0KLWGcA0xz9wHAtGT7U8ysB3AesDcwCDivUWL5lbvvAOwO7Gdmwxu99GZ33y25jSswzqLYe2/YZBNVS4lIdSo0YYwExiePxwNH5DjmEGCquy9x96XAVOBQd3/f3f8O4O4rgaeALQuMp6TWWw8OOyzmlVq9Ou1oRETKq9CE0dvdFwEk95vmOGYL4NVG2wuTff9hZt2ALxKllAZHmdksM5tgZn0LjLNosllYsgQeeyztSEREyqvZhGFm95vZszluI/O8huXY95/hb2ZWA9wIXObuLyW77wK2dvfPAfezphSTK75TzKzOzOrq6+vzDKn1hg2LkoaqpUSk2jSbMNz9YHffJcdtIrDYzPoAJPdv5jjFQqBxCWFL4PVG21cB89z90kbXfNvdG6b6uxrYcx3xXeXute5e26tXr+Z+nIJ16waZjGavFZHqU2iV1CRgVPJ4FDAxxzH3AsPMrHvS2D0s2YeZ/QLYGPhB4xc0JKHE4cDcAuMsqmwWZs+GV15JOxIRkfIpNGFcCAw1s3nA0GQbM6s1s3EA7r4EOB94Mrn93N2XmNmWwLnATsBTa3WfHZN0tZ0JjAFOKjDOompYVEmlDBGpJubtaDa92tpar6urK/l13GG77WDAAK3EJyJtn5nNcPfa5o7TSO9WMItR3w88AO+9l3Y0IiLloYTRStlsrMD3wANpRyIiUh5KGK10wAHQtau614pI9VDCaKVOneCQQ7SokohUDyWMAmSz8Prr8MwzaUciIlJ6ShgFGJ5MlahqKRGpBkoYBejdGwYN0ngMEakOShgFymbhiSdg8eK0IxERKS0ljAJls9HofffdaUciIlJaShgF2m032HxztWOISPunhFGghlHf990HK1emHY2ISOkoYRRBNgsrVsD06WlHIiJSOjVpB9AeDBkC668f1VJDhqQdjUj79NZbMHPmmtu228JPf5p2VNVFs9UWyWGHwQsvwLx5UU0lIq2zalX8HTVODjNnxiDZBhtsAB9/DMuWxRQ9Uph8Z6tVCaNIRoyInlIvvADbb592NCJtw7JlaxLCrFlx/+yz8OGH8XxNDey0U5TcP/c5GDgwbk8/DYceCo8+CkOHpvszVBMljCIZMQK+970YxKeEURznnQfdu8OYMdBBrW1t2iefwIsv/nep4d//XnNMz56RDL7znTWJYccdY962tQ0eHJ+J6dOVMMqp4CopM+sB3AxsDbwMHOPuS3McNwr4cbL5C3cfn+z/B9AH+CB5bpi7v2lm6wPXE+t5vw18xd1fXlcsaVZJAey6K/TqpSnPi+Gpp2DPZCX3bBbGj4cePdKNSfKzYkUsYdw4McyevWbtmA4d4p+qhqTQcOvTp2XVuXvuCZ/5DPz976X5OapJOaukzgGmufuFZnZOsn32WsH0AM4DagEHZpjZpEaJ5Xh3X/ub/mRgqbt/1syOBS4CvlKEeEsmm4Vf/SqK2d26pR1N2zZ2LGy0EfzkJ3DuubDHHnDrrbDXXmlHJg3cY137tUsNL7645phu3SIZnHzymiqlnXeGzp0Lv34mA3/8Y3Rnz1UKkRJw94JuwPNAn+RxH+D5HMccB/yx0fYfgeOSx/8AanO85l5gcPK4BniLpETU1G3PPff0ND30kDu433xzqmG0eS++6N6hg/uZZ8b244+7b7WVe8eO7pdf7v7JJ+nGV83+/W/3MWPcMxn3jTeOzzu4m7kPGOB+9NHu55/vPmmS+yuvlPZ3NWFCXPuRR0p3jWoB1Hke3/fFKGH0dvdFSfJZZGab5jhmC+DVRtsLk30N/mRmq4HbiOoqb/wad19lZsuBTZLE8R9mdgpwCsBWW21VhB+n9fbZJ6pNJk+GY45JNZQ27de/jsbOH/wgtgcNikbOE0+E//1feOghuPrqKIFI+bzzTszQPH9+VAd99atrqpN22aX8vZX23z/up0+PNg0pvbwShpndD2yW46lz87xOrprJhsaT4939NTPbiEgYJxBtF+t6zZod7lcBV0G0YeQZT0mst178Qd19N6xeHdvSMm++CddeCyecEFOuNOjRAyZNiqqqc8+NNUgmTIgvKim91avh+OPhuefg3nsrY7xR796w3XaRMM46K+1oqkNefU/c/WB33yXHbSKw2Mz6ACT3b+Y4xUKgb6PtLYHXk3O/ltyvAP4KDFr7NWZWA2wMLGnpD1hu2WwMMHriibQjaZuuuCLWSj/jjP9+rkMHOOccmDYt2okGDYrGcCm9c8+NkvNvf1sZyaLBAQfAww9HLywpvWJ0VpwEjEoejwIm5jjmXmCYmXU3s+7AMOBeM6sxs54AZtYRyALP5jjv0cADSVVVRTvkkChZaDLClnv33UgYI0fCDjs0fdyBB0YJY++94aSTYPRo+OCDpo+XwvzlL3DRRfCtb0WX10qSycDSpTBnTtqRVIdiJIwLgaFmNg8YmmxjZrVmNg7A3ZcA5wNPJrefJ/vWJxLHLOAZ4DXg6uS81wCbmNl84DSi91XF69496laVMFpu3Lj44z/77OaP3WwzmDoVfvQjuOaaqMOeN6/0MVabxx+PhHzggXDZZZU3i0EmE/eax608NDVICVx8cdSpvvIKpNwO32Z8/HHMDdS/P/zzny177ZQp0ebx8cfwpz/BUUeVJsZqs3BhdGPu3DmqWHv2TDui/+YOfftG4rjxxrSjabvyHYeh8bMlkM3G/ZQp6cbRltx0E7z6ausaLw87LAb67bgjHH00nHqqppov1PvvwxFHRDXhXXdVZrKAKPFkMlHCaEf/+1YsJYwS2GEH2GYbVUvlyz16P+2yS3z5t0a/fvGl8b//C5deGlUor77a7MskB3f4+tcjCd94Ywy0q2SZDLz2Grz8ctqRtH9KGCVgFqWMadPiPzVZtylTYsK5s84qrI68U6eoZ7/55piKYvfd4Z57ihdntbjggngPL7xwTWm5kqkdo3yUMEokm40ZNzWvVPPGjo166GOPLc75jjkGZsyIcRyHHRZrJqxeXZxzt3d33BHTsZxwApx5ZtrR5GfnnaOziRJG6SlhlMgBB8CGG8bstdK0xx6DBx+E006Djh2Ld97ttotzjxoF558f3Z0XLy7e+dujmTPhf/4nuitfdVXl9YhqSocOsN9+ShjloIRRIuuvD8OGRTuGGuOadtFF8d/h6NHFP3eXLtFr6pprYnDX7rvrS6Upb74Jhx8ev4s77ogFitqSTAaefz5+DikdJYwSymaja+KsWWlHUpmeew4mTox1REo5D9HXvx6ljQ03hC98IarAlMTX+Ogj+NKXoL4+fh99+qQdUcs1tGM89FC6cbR3Shgl1NDjR72lcvvVr6Ik9r3vlf5aAwdGu8aRR8bAwCOOiEGC1c49Rm8//HCUxhrWIGlr9twzxouoBFlaShgltNlmMfBJCeO/vf46/PnP8d//prnmNy6Bz3wGbrkl5kOaMiXW2JgxozzXrlS//W1M9viTn8BXKnq1mXXr1CnaXpQwSksJo8RGjIjpFerr046kslx6KaxaBaefXt7rmsWSr9OnR8+pffeF3/++Oquo7r033v8jj4Sf/SztaAqXycQ0+CtWpB1J+6WEUWLZbHwZ3X132pFUjuXL4Q9/iO6v22yTTgz77BMD0w46KKpkjj8+RjVXi+eeixLFrrvC9de3jzXTM5mYtfbRR9OOpP1qBx+Tyrb77tGIqGqpNf7wh/gvMO01DHr2jG7Pv/hFDFTba6/qmPV06dLoEdWpUzRyl3vho1IZPDhmila1VOkoYZRYhw5RLXXvvZrfCGIw46WXwtChkUzT1qFDrPUwdSosWRJrbPzlL2lHVTqrVkXJ7uWXo/tsv35pR1Q8Xbuq63SpKWGUQTYby1uqy180dL/xRn5TmJfTQQdF/feee8Yo529+M5Jbe3P66XD//VHK22+/tKMpvkwm2gw/+ijtSNonJYwyGDIkiv/VPup79eroSrvHHvEFXWk23zymcjn77BjpvO++8OKLaUdVPFdfHXNtnXpq9E5rjzKZSPTV3vutVApKGGbWw8ymmtm85L57E8eNSo6ZZ2ajkn0bmdkzjW5vmdmlyXMnmVl9o+dKMA64fLp2jQFj1d6OMXEivPBCfCFX6rQTNTUx6d6kSbBgQZQ47rgj7agK9+CD0bh/6KExcLG92n//uFe1VGkUWsI4B5jm7gOAaeRYFc/MegDnAXsT63WfZ2bd3X2Fu+/WcANeAW5v9NKbGz0/rsA4U5fNxpflCy+kHUk63GMakG23bRsLHH3xi9GL6rOfjVHQJ5wQC2K1RQsWxHu+7bYxXXlNTdoRlU6vXrG8gBJGaRSaMEYC45PH44EjchxzCDDV3Ze4+1JgKnBo4wPMbACwKdBuf80jRsR9tVZL/fOfsWrbGWdET5a2oH//GAF9zjlw660xoeHpp0fjeFuxYkX0iFq1KhZC6tYt7YhKL5OJ39snn6QdSftTaMLo7e6LAJL7XGN2twAaL2WzMNnX2HFEiaLx8KmjzGyWmU0ws74Fxpm6/v1jGuZqrZYaOzZGdI8alXYkLbP++vDLX8Z64V/9KlxySfynPnYsfPBB2tGt2yefxOyzc+fGCPcBA9KOqDwyGVi2LNZYkeJqNmGY2f1m9myO28g8r5GrtnrtcbXHAo1X5L0L2NrdPwfcz5pSTK74TjGzOjOrq6/w4dQjRkRd8jvvpB1Jec2aFQMXx4yJ+X7aor59Y66lmTOjMfzss2H77eG66yp3rY0f/zjaYi65JLoxVwstqFQ6zSYMdz/Y3XfJcZsILDazPgDJfa7JhRcCjUsIWwKvN2yY2UCgxt3/06/B3d9294aOcVcDTU6J5u5XuXutu9f26tWruR8nVdlsVA3cd1/akZTX2LHR8P+d76QdSeF23TWqFR94AHr3hq99Lfr+3313ZU0v8te/RsnoG98oz+SOlaRfP9hySyWMUii0SmoS0FDJMAqYmOOYe4FhZtY96UU1LNnX4Dg+XbpoSD4NDgfmFhhnRRg8ONYbqKZqqVdegZtuglNOiZ+9vfjCF6K//003wXvvxczEQ4ZAXV3akUVb0de/Hot4XXFF5fZIKxWzKGVMn15ZSbw9KDRhXAgMNbN5wNBkGzOrNbNxAO6+BDgfeDK5/TzZ1+AY1koYwBgzm2NmM4ExwEkFxlkRampg+PCYKbVaGuR+85v4A/7BD9KOpPg6dIj5mObOjVlfZ8+O6UWOOw5eeimdmF57LaZu79MHbrstxv9Uo0wmZkResCDtSNoZd283tz333NMr3Q03uIP7Y4+lHUnpvfWWe5cu7qNGpR1JeSxf7n7uue6dO7t37Og+Zox7fX35rv/+++61te5du7rPmlW+61ai2bPj7+y669KOpG0A6jyP71iN9C6zQw+N/0yroVrqiivg/ffhzDPTjqQ8PvOZmMhw/nw46aT4+bfdFv7f/4v3oZTcoxpqxgy44YZoa6lmO+0UVaBqxyguJYwy69Ej5vBp7wnj/ffh8sujoX/nndOOprw23zymFpk9Gw48MCY3HDAg1hZftao01/zlL6M95YILYtxFtevQIUZ9K2EUlxJGCrJZeOaZWO+7vbr2Wnj77cqbZLCcdtoppkN58EHYaisYPTqWir3rruI2xk6cGEnpq1+NQYYSMpmYWWHx4rQjaT+UMFKQzcZ9ex31vWoV/PrXMV6hYW6fapbJwCOPwIQJ8PHHUQL4/Oejl1WhZs2KxZ/22gvGjau+HlHr0jAeQ7NEF48SRgp23BG23rr9Joxbb431Fqq5dLE2s5jPac4cuPJKeP75WPXvy1+OUeStUV8fyWfjjeHOO9vuoMhS2WOPeE9ULVU8ShgpMItSxv33V/70Ei3VMMngjjuuKUnJGh07wre/HQ3j550XA/522ikG172Za9hrE1aujAS0eHEki803L13MbVWnTpGUlTCKRwkjJdlsJIu//z3tSIrrvvti+owzz2wf60SXykYbwc9+FonjG9+IBY223RZ+/vPm1xZ3h+9+N74Ir702qqMkt0wm2gurbTqeUtGfdEo+/3nYcMP211vqooviv93jj087krZhs82iimrOHBg2LEodn/1sJJCPP879mssvj/aKc8+NQYLStEwmBsk++mjakbQPShgp2WCDmBBu8uT2M33Bk09GienUU6t3hHFrbb99jMx+5JHogvvtb8dYijvu+PTn47774v094ogojci67bNPTKevaqniUMJI0YgR8Oqr7Wca5rFjowH2lFPSjqTtGjw4uuFOnBhtXV/6UvQ0e/jh6CL6la/EuJY//1lVfvno2jUav5UwikMfuRQddljct4dqqXnz4j/k73wnRjxL65lF76fZs2MA4IIFkTQGDYr5yCZNii9CyU8mE12YP/qo+WNl3ZQwUrT55rFmdHtIGL/6VVRDjRmTdiTtR01NNIjPmxdTjmy2Gdx+e3TJlvxlMpEsKmEm4bZOCSNl2Ww0yL31VtqRtN4bb8D48bGa3mabpR1N+7PhhtHA/dxzawajSf4aBo+qWqpwShgpy2ajUfPuu9OOpPUuuyzGBZxxRtqRiPy3nj1jXJASRuGUMFK2xx6xcltbHfX9zjvRLfSoo6pnzWhpezKZ6DhQqcvpthUFJwwz62FmU81sXnKfc101M7vHzJaZ2eS19vc3s8eT199sZp2S/esn2/OT57cuNNZK1KFD9Ja6556m+91XsquvhuXL4ayz0o5EpGmZTHxO20uPxLQUo4RxDjDN3QcA05LtXC4GTsix/yLgkuT1S4GTk/0nA0vd/bPAJclx7VI2Gx/mhx9OO5KWWbkSLrkklivVaGOpZA1tP6qWKkwxEsZIYHzyeDxwRK6D3H0asKLxPjMz4CBgQo7XNz7vBGBIcny7c/DB0cOorfWWuuGGWBJUkwxKpevXD/r2VcIoVDESRm93XwSQ3G/agtduAixz94ZlZRYCWySPtwBeTc67ClieHN/ubLRRLLTTltoxPvkkBuoNHBhTWohUukwmEkZ7mVkhDXklDDO738yezXEbWeD1c5UYPI/nGsd2ipnVmVldfX19geGkZ8SI6DY5f37akeRn8uSI96yztAaDtA2ZDCxaBC+9lHYkbVdeCcPdD3b3XXLcJgKLzawPQHLfgkmaeQvoZmY1yfaWwOvJ44VA3+S8NcDGwJIcsV3l7rXuXturV68WXLqyjBgR922llHHRRTGA7Jhj0o5EJD9qxyhcMaqkJgGjksejgIn5vtDdHfg7cHSO1zc+79HAA8nx7dK220Zf8bvuSjuS5j30UEySd/rpMRpZpC3YcUfo0UMJoxDFSBgXAkPNbB4wNNnGzGrNbFzDQWY2HbiVaLxeaGaHJE+dDZxmZvOJNoprkv3XAJsk+0+j6d5X7cZRR8G0adEmUIzlO0tl7FjYZBP4+tfTjkQkfx06xKhvJYzWK/j/Q3d/GxiSY38dMLrRds5JDdz9JWBQjv0fAl8uNL625Mc/hu7d4Ze/jGmZv/jFmMJ6t93SjmyNOXOiFPSzn0GXLmlHI9IymUxM3vjGG5rGpjU00ruCrL8+nHZaNMpdcEH8J7T77tFOMHdu2tGFiy+ORPG976UdiUjLNbRjPPRQunG0VUoYFWijjeBHP4pprX/845hnapdd4MQT4cUX04vr1Vdj7MXo0VElJdLW7LFH/MOjaqnWUcKoYN26wfnnR4njtNPg1lthhx3gm9+ML+9yu/TS6MN+2mnlv7ZIMXTsGNW9Shito4TRBvTqFVVBL70E3/oW/OlPse7z978fdbHlsHRpLOZz7LExalakrcpkYObMmDhTWkYJow3p0wcuvzwW1DnxRPjd72CbbWJqjrffLu21r7wS3n1XkwxK25fJxEwFjzySdiRtjxJGG9SvX8wSO3durPl88cXQvz+cd15MYlhsH3wAv/0tDB8On/tc8c8vUk777BPjh1Qt1XJKGG3YgAHwl7/E2s/DhkUX3P794cIL4b33ined8eOhvl6lC2kfNtwwGr+VMFpOCaMd2HlnmDABZsyAffeFH/4wqqouvRQ+/LCwc69eHet1DxoEn/98ceIVSVsmE4NjC/37qDZKGO3IHnvEpICPPAK77gqnnhqN43/4Q6xd0Rq33RZdec8+W5MMSvuRycTfxJNPph1J26KE0Q4NHgz33w8PPBDtHd/+Nmy/PVx3Haxa1ezL/8M9pgHZbjsYWei8xCIVZP/9417VUi2jhNGOfeELMaJ1ypSYdO1rX4sBgDffHL1EmvPAA1HNdcYZsN56pY9XpFw22QR22kkJo6WUMNo5s+jdVFcHt98evUOOPTbmp5o4cd2LyVx0Ucy3c0KuhXVF2rhMJqpvV69OO5K2QwmjSpjBkUfGgKW//jUa+444Ihqz7733vxPHU0/B1Knwgx/ABhukE7NIKWUyMXhv1qy0I2k7lDCqzHrrwXHHwb/+BddeG91lDz00ekA9+OCa4y6+OOa0+uY304tVpJS0oFLLKWFUqZqaaNN44YUYxf3ii5E0hg2LOatuuSWmIenWLe1IRUpjq63ipoSRPyWMKtepU/Simj8ffv1rePrpmE69piaqo0Tas0wmEkb7XcuzuApKGGbWw8ymmtm85L57E8fdY2bLzGzyWvtvMLPnzexZM7vWzDom+w80s+Vm9kxy+2khcUrzOneOWWgXLIjG7iuvhM03TzsqkdLKZGDx4viHSZpXaAnjHGCauw8AptH0MqoXA7n62twA7ADsCnSm0Qp9wHR33y25/bzAOCVPXbvGFCAnn5x2JCKlp3aMlik0YYwExiePxwNH5DrI3acBK3Lsn+IJ4AlgywLjERHJ2447xpgMJYz8FJowerv7IoDkftPWnCSpijoBuKfR7sFmNtPM7jazndfx2lPMrM7M6urr61tzeRGpUmYx6lsJIz/NJgwzuz9pY1j7VszJIq4EHnT3hl/bU0A/dx8IXA7c2dQL3f0qd69199pevXoVMSQRqQaZTPQSXLQo7UgqX01zB7j7wU09Z2aLzayPuy8ysz7Amy0NwMzOA3oB/+nx7+7vNHo8xcyuNLOe7v5WS88vIrIujdsxjjkm3VgqXaFVUpOAUcnjUcDElrzYzEYDhwDHufsnjfZvZhZzo5rZoCTOEq8pJyLVaPfdoUsXVUvlo9CEcSEw1MzmAUOTbcys1szGNRxkZtOBW4EhZrbQzA5JnvoD0Bt4dK3us0cDz5rZTOAy4NikYVxEpKg6dowZnpUwmtdsldS6uPvbwJAc++to1EXW3TNNvD7n9d39CuCKQmITEclXJgP/93+wbJlmN1gXjfQWkaqXycRo70ceSTuSyqaEISJVb599YjocVUutmxKGiFS9Ll1gzz2VMJqjhCEiQlRLPflkrBUjuSlhiIgQCWPlSnjiibQjqVxKGCIiwH77xb2/wnfqAAAL5UlEQVSqpZqmhCEiQkxCuPPOShjrooQhIpLIZKJr7erVaUdSmZQwREQSmQysWAEzZ6YdSWVSwhARSWhBpXVTwhARSfTtC/36KWE0RQlDRKSRAw6IhKHpTv+bEoaISCOZDLz5Jsybl3YklUcJQ0SkEbVjNE0JQ0Skke23h169lDByKShhmFkPM5tqZvOS++5NHHePmS0zs8lr7b/OzBYkiyc9Y2a7JfvNzC4zs/lmNsvM9igkThGRfJnB/vsrYeRSaAnjHGCauw8ApiXbuVwMnNDEc2e6+27J7Zlk33BgQHI7Bfh9gXGKiOQtk4GXXoLXX087kspSaMIYCYxPHo8Hjsh1kLtPA1a08LzXe3gM6GZmfQqKVEQkT22tHaNcPboKTRi93X0RQHK/aSvOcUFS7XSJma2f7NsCeLXRMQuTfSIiJbfbbtC1a9tIGB9/DEcdBbffXvprNZswzOx+M3s2x21kEa7/Q2AHYC+gB3B2w2VzHJszh5rZKWZWZ2Z19fX1RQhJRKpdTQ0MHlz5CcMdvvENuOMOWLKk9NdrNmG4+8HuvkuO20RgcUNVUXL/Zksu7u6Lkmqnj4A/AYOSpxYCfRsduiWQszbR3a9y91p3r+3Vq1dLLi8i0qRMBmbPhmXL0o6kaT/5CYwfDz/7GYweXfrrFVolNQkYlTweBUxsyYsbJRsj2j+ebXTeE5PeUvsAyxuqvkREyiGTif/gH3447Uhy+/3v4YILooTx05+W55qFJowLgaFmNg8YmmxjZrVmNq7hIDObDtwKDDGzhWZ2SPLUDWY2G5gN9AR+keyfArwEzAeuBr5TYJwiIi2y997QsWNlVkvdcQd897vwxS/ClVdGV+ByMG9HE6bU1tZ6XV1d2mGISDux777xZVxJpYyHHoKDD4bdd4dp06BLl8LPaWYz3L22ueM00ltEpAmZDDz5JHzwQdqRhH/9Cw4/PGbUveuu4iSLllDCEBFpQiYT3VafeCLtSOC11+DQQ2H99eGee6Bnz/LHoIQhItKE/faLKqm02zGWL4fhw6PH1pQp0L9/OnHUpHNZEZHK17077LJLugnjo4/giCNg7ly4++5ou0iLShgiIuuQycAjj8CqVeW/9iefwIknwj/+AdddF43daVLCEBFZh0wG3n0XZs4s/7XPOANuuQXGjoXjjy//9demhCEisg5pTUT461/DJZfA978fiaMSKGGIiKzDFltEI3M5E8Zf/xpJ4stfht/8pnwD85qjhCEi0oxMJhJGOcY5T5sGJ50En/88XH89dKigb+kKCkVEpDJlMlBfDy+8UNrrPPMMHHlkLBN7552wwQalvV5LKWGIiDSjHO0YL78cYy023ji6z3brVrprtZYShohIM7bbDjbdtHQJ4+23YxT3hx/GKO4ttyzNdQqlgXsiIs0wg/33L03CeP/9mHX25Zdh6lTYeefiX6NYVMIQEclDJgMLFsScTsWyahUcdxw89lj0jGqo+qpUShgiInkodjuGe6xpMWkSXH45fOlLxTlvKRWUMMysh5lNNbN5yX33Jo67x8yWmdnktfZPN7NnktvrZnZnsv9AM1ve6LkyrSclIpLbwIHQtWvxEsYvfgFXXQU//GEkjrag0BLGOcA0dx8ATEu2c7kYOGHtne6ecffd3H034FHg9kZPT294zt1/XmCcIiIFqamJBZWKkTCuuSaWVT3xxFhmta0oNGGMBMYnj8cT63L/F3efBqxo6iRmthFwEHBngfGIiJRMJgPPPgtLl7b+HH/7G3zzm3DIITBuXOWM4s5HoQmjt7svAkjuN23leY4kSirvNNo32MxmmtndZlbB/QZEpFpkMtH20NolW594Ao45BnbbDSZMiDXD25Jmu9Wa2f3AZjmeOreIcRwHjGu0/RTQz93fNbPDiJLHgCbiOwU4BWCrrbYqYkgiIp82aFB8yU+fDtlsy177wgswYgRstlmUMrp2LU2MpdRswnD3JmdgN7PFZtbH3ReZWR/gzZYGYGabAIOIUkbDNd9p9HiKmV1pZj3d/a0c8V0FXAVQW1tbhpleRKRade4Me+3V8naMxYtjYB7EwLzevYsfWzkUWiU1CRiVPB4FTGzFOb4MTHb3Dxt2mNlmZlGzZ2aDkjjfLjBWEZGCZTJQVwcffJDf8StWwGGHRdL4299gQM66krah0IRxITDUzOYBQ5NtzKzWzP5TxWRm04FbgSFmttDMDml0jmOBG9c679HAs2Y2E7gMONa9HPNEioisWyYDH38Mjz/e/LErV8LRR8fiS7feGlVabVlBU4O4+9vAkBz764DRjbabHL/o7gfm2HcFcEUhsYmIlMJ++0XPpunT4cADmz7OHUaPhvvug2uvjVJGW6eR3iIiLdCtG+y6a/PtGD/6Efz5z3D++fC1r5UntlJTwhARaaFMBh59NOaCyuWKK+DCC2O8xbnF7E+aMiUMEZEWymTg3XdjwaO13XYbjBkDI0fC737XtgbmNUcJQ0SkhZqaiHD6dDj+eBg8GG68EdZbr/yxlZIShohIC22+OWyzzacTxpw5cPjh0L9/zEDbuXN68ZWKEoaISCtkMvDQQ9EbauHCGJjXuXMMzNtkk7SjKw0lDBGRVshkoL4+xmMMHw7Ll8da3P36pR1Z6WiJVhGRVmhoxxg+HN57L0oWAwemG1OpqYQhItIKAwbAppvCsmVw/fVw0EFpR1R6KmGIiLSCGYwdCx06wLHHph1NeShhiIi00qhRzR/TnqhKSkRE8qKEISIieVHCEBGRvChhiIhIXpQwREQkL0oYIiKSFyUMERHJixKGiIjkxdw97RiKxszqgVfSjqNAPYG30g6iguj9+DS9H2vovfi0Qt6Pfu7eq7mD2lXCaA/MrM7da9OOo1Lo/fg0vR9r6L34tHK8H6qSEhGRvChhiIhIXpQwKs9VaQdQYfR+fJrejzX0Xnxayd8PtWGIiEheVMIQEZG8KGGkyMz6mtnfzWyumc0xs+8n+3uY2VQzm5fcd0871nIxs/XM7Gkzm5xs9zezx5P34mYz65R2jOViZt3MbIKZPZd8RgZX+Wfj1OTv5Fkzu9HMNqiWz4eZXWtmb5rZs4325fwsWLjMzOab2Swz26NYcShhpGsVcLq77wjsA3zXzHYCzgGmufsAYFqyXS2+D8xttH0RcEnyXiwFTk4lqnT8FrjH3XcABhLvS1V+NsxsC2AMUOvuuwDrAcdSPZ+P64BD19rX1GdhODAguZ0C/L5YQShhpMjdF7n7U8njFcQXwhbASGB8cth44Ih0IiwvM9sSGAGMS7YNOAiYkBxSTe/FZ4ADgGsA3H2luy+jSj8biRqgs5nVAF2ARVTJ58PdHwSWrLW7qc/CSOB6D48B3cysTzHiUMKoEGa2NbA78DjQ290XQSQVYNP0IiurS4GzgE+S7U2AZe6+KtleSCTUarANUA/8KamiG2dmG1Klnw13fw34FfBvIlEsB2ZQvZ8PaPqzsAXwaqPjiva+KGFUADPrCtwG/MDd30k7njSYWRZ4091nNN6d49Bq6dZXA+wB/N7ddwfeo0qqn3JJ6udHAv2BzYENiaqXtVXL52NdSvZ3o4SRMjPrSCSLG9z99mT34oYiZHL/ZlrxldF+wOFm9jJwE1HVcClRnK5JjtkSeD2d8MpuIbDQ3R9PticQCaQaPxsABwML3L3e3T8Gbgf2pXo/H9D0Z2Eh0LfRcUV7X5QwUpTU0V8DzHX33zR6ahIwKnk8CphY7tjKzd1/6O5buvvWRGPmA+5+PPB34OjksKp4LwDc/Q3gVTPbPtk1BPgXVfjZSPwb2MfMuiR/Nw3vR1V+PhJNfRYmAScmvaX2AZY3VF0VSgP3UmRm+wPTgdmsqbf/EdGOcQuwFfGH8mV3X7vBq90yswOBM9w9a2bbECWOHsDTwP+4+0dpxlcuZrYb0QGgE/AS8DXin7yq/GyY2f8BXyF6Fz4NjCbq5tv958PMbgQOJGakXQycB9xJjs9CklCvIHpVvQ98zd3rihKHEoaIiORDVVIiIpIXJQwREcmLEoaIiORFCUNERPKihCEiInlRwhARkbwoYYiISF6UMEREJC//HzbqVoeTRIr2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x3300528128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 绘制不同聚类数目的模型性能， 找到最佳模型\n",
    "plt.plot(Ks, np.array(CH_scores), 'b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最佳的K值是10，在10附近再试试看是不是有更好的值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K-means begin with 5 clusters\n",
      "CH_score: 0.18870415532631862\n",
      "K-means begin with 6 clusters\n",
      "CH_score: 0.11313165718804834\n",
      "K-means begin with 7 clusters\n",
      "CH_score: 0.17582010267437606\n",
      "K-means begin with 8 clusters\n",
      "CH_score: 0.042239916967105585\n",
      "K-means begin with 9 clusters\n",
      "CH_score: 0.06368888437045797\n",
      "K-means begin with 10 clusters\n",
      "CH_score: 0.04116859444177382\n",
      "K-means begin with 11 clusters\n",
      "CH_score: 0.10934080630736519\n",
      "K-means begin with 12 clusters\n",
      "CH_score: -0.07955099588307381\n",
      "K-means begin with 13 clusters\n",
      "CH_score: 0.10510935007582564\n",
      "K-means begin with 14 clusters\n",
      "CH_score: -0.013461979410547633\n",
      "K-means begin with 15 clusters\n",
      "CH_score: 0.014042538867224116\n"
     ]
    }
   ],
   "source": [
    "# 调整超参数（聚类数目K）的搜索范围\n",
    "CH_scores = []\n",
    "Ks = [5,6,7,8,9,10,11,12,13,14,15]\n",
    "for K in Ks:\n",
    "    ch = K_cluster_analysis(K, eventContMatrix.toarray())\n",
    "    CH_scores.append(ch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x33009a6be0>]"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XuclHX5//HXBctJxBNCKmAcBfGQ6AJjKpaIoqmDicpiX/XXwbLMMjM1zWOapqXmIbWjmabmEY9IWloxAQsqgiALSBxERPGIymH5/P64ZnJ22V12d+6Ze2b2/Xw89jG7s/fc97Ue5j3352ghBERERDLaxV2AiIgUFwWDiIjUoWAQEZE6FAwiIlKHgkFEROpQMIiISB0KBhERqSOSYDCzsWb2qpktNLPzGvj9D8zsFTObbWbPmNlns353ipnVpL9OiaIeERFpPct1gpuZtQcWAGOA5cAMoCqE8ErWMV8EpoUQPjKz04EvhBBONLMdgGqgEgjATGC/EMI7ORUlIiKtVhHBOUYAC0MIiwHM7B4gCfwvGEIIf886/j/AV9LfHw5MCSGsSb92CjAW+EtTF9xxxx1D3759IyhdRKTtmDlz5lshhB5bOi6KYOgFLMv6eTkwsonjvwY82cRre23pgn379qW6urqFZYqItG1m9t/mHBdFMFgDzzXYPmVmX8GbjQ5uxWtPA04D2HXXXVtepYiINEsUnc/LgT5ZP/cGXq9/kJkdClwAHBNCWNeS1wKEEG4PIVSGECp79NjinZCIiLRSFMEwAxhkZv3MrCMwAZiUfYCZDQNuw0PhzaxfTQYOM7PtzWx74LD0cyIiEpOcm5JCCBvN7Az8Db098PsQwlwzuwyoDiFMAq4Btgb+amYAS0MIx4QQ1pjZ5Xi4AFyW6YgWEZF45DxcNQ6VlZVBnc8iIi1jZjNDCJVbOk4zn0VEpA4Fg4iI1NFmgmHTJvjd7+DBB+OuRESkuEUxj6Fk/PrX8OabcOSR0Llz3NWIiBSnNnPH0K4dXHMNLFsGN94YdzUiIsWrzQQDwBe/6HcLV1wBb78ddzUiIsWpTQUDwNVXwwcfwJVXxl2JiEhxanPBsOeecOqpcNNNsGRJ3NWIiBSfNhcMAJdeCu3bw4UXxl2JiEjxaZPB0Ls3nHUW3HUXzJoVdzUiIsWlTQYDwI9+BN27wznnQAmuCiIikjdtNhi23RYuugiefRYmaz1XEZH/abPBAPCtb0H//n73UFsbdzXRe/dd2H13+NWv4q5EREpJmw6Gjh3hZz+Dl1+GO++Mu5ronX02zJ8Pv/993JWISClp08EAcPzxMHy4j1D6+OO4q4nO5MkeCP36wUsvaWiuiDRfmw8GM18qY8UKuOGGuKuJxvvvwze+AUOGwKT0XnqTJjX9GhGRjDYfDAAHHwxHH+3NSm+9FXc1uTv3XFi+3O8Y9tzT+xkeeSTuqkSkVCgY0q66Cj78EH7607gryc2zz8Ktt/o8jf339+eSSXjuOXjnnXhrE5HSoGBIGzoUvvY1uOUWWLQo7mpa58MP4etfh4ED4fLLP31+3DgfdfX44/HVJiKlQ8GQ5ZJLoEMHuOCCuCtpnQsugNde8w2Jttrq0+eHD4edd1Zzkog0j4Ihyy67+BDPe++FGTPirqZl/vUv32fijDNg1Ki6v2vXzvtQnnoK1q2Lpz4RKR0KhnrOOQd69CitpTI++gi++lXo29c70BuSTHpT07PPFrQ0ESlBCoZ6unWDiy/2ztonnoi7mua5+GKoqYHf/Aa23rrhYw45BLp2VXOSiGyZgqEBp50Ggwb5UhkbN8ZdTdOmTYNf/tJrHj268eM6d4axY30+w6ZNhatPREqPgqEBHTp4k8wrr8Add8RdTeM++QT+3/+DXr18kt6WjBsHK1eWXv+JiBSWgqERX/4yJBK+AuvatXFX07DLL4d58+D222GbbbZ8/JFH+gZFak4SkaYoGBqRWSrj9dfh+uvjrmZzM2f6/tWnnupNRM2xww4+YknBICJNUTA04cADvfnl6qvhzTfjruZT69d7E1LPnt6/0BLJpDeRLVyYn9pEpPRFEgxmNtbMXjWzhWZ2XgO/H2Vms8xso5mNr/e7WjN7Mf1VdEu9/exnPhw0eyZx3DJLhd96K2y/fctem0z6o+4aRKQxOQeDmbUHbgaOAIYCVWY2tN5hS4FTgbsbOMXHIYR90l/H5FpP1IYM8ZVKb73Vh4TGbfZsX89p4kQ4phX/tPr2hb33VjCISOOiuGMYASwMISwOIawH7gGS2QeEEJaEEGYDJTlQ8uKLoVMn+PGP461jwwZvQtphh9x2ZRs3Dv79b1i9OrraRKR8RBEMvYBlWT8vTz/XXJ3NrNrM/mNm4yKoJ3I77eQzoe+/H/7zn/jquPZamDXLF/rr3r3150kmfS7DY49FV5uIlI8ogsEaeK4li0nsGkKoBCYC15vZgAYvYnZaOkCqV8fwUffss+Ezn4lvqYxXXvFF/saPh+OOy+1cw4ZBnz5qThKRhkURDMuBPlk/9wZeb+6LQwivpx8XA/8AhjVy3O0hhMoQQmWPHj1aX20rbb01XHqpL1ZX6N3Qamt9LaRu3eDmm3M/n5n3Tzz9tHesi4hkiyIYZgCDzKyfmXUEJgDNeus0s+3NrFP6+x2BA4BXIqgpL772NRg8GM47r7BLZVx/vS99ceONPkQ1Csmk73H9t79Fcz4RKR85B0MIYSNwBjAZmAfcF0KYa2aXmdkxAGY23MyWA8cDt5nZ3PTLdweqzewl4O/AVSGEog2Gigrf6W3+fN82sxAWLIALL/RP+BMmRHfegw/22dJqThKR+iyUytrSWSorK0N1dXUs1w4BDjrId3mrqWl8NdMobNrkb+Bz5ngfw847R3v+iRP9jmHlSl8qQ0TKm5nNTPfpNkkzn1sos1TGG2+0fNZxS918s/dpXH999KEA3py0ejWkUtGfW0RKl4KhFfbf30cG/fznsGpVfq6xeLH3ZRxxBJx8cn6uccQRvpKsmpNEJJuCoZWuvNK3ybz00ujPvWkTfP3r3rxz221+l5IP22wDX/yiB0MJtiiKSJ4oGFppt93gm9/0Ja9ffTXac//mN/D3v8MvfuHzDfIpmfS+kvnz83sdESkdCoYcXHQRdOkC558f3TmXLvVJdKNH+11DvmXWW1JzkohkKBhy0LMnnHsuPPSQrz2UqxB8wb5Nm+C3v81fE1K23r1hv/0UDCLyKQVDjs46y0cMRbFUxh//6LORr77aV0EtlHHjfA2olSsLd00RKV4Khhx17QqXXeZDPh96qPXnWbHCQ2bUKDj99Ojqa47MHg2PPlrY64pIcVIwRODUU2HoUB9eumFDy18fAnzrW74z2+9+B+0K/G9lzz2hXz81J4mIUzBEoKLCm39qanxEUUvdfbcvgX3FFTBwYPT1bYmZ3zU88wx8+GHhry8ixUXBEJEvfcmbgS65BD74oPmve+MNOPNMnzR35pl5K2+LkkmflzF5cnw1iEhxUDBEJLNUxurVvqFOc4QA3/kOrF3ri/LFuV7RgQf6znBqThIRBUOERoyAE07wYGjOCJ/774cHH/TZ00OG5L++plRUwFFHeZNWa/pJRKR8KBgiduWV/sZ6ySVNH7d6td8tVFb67nDFIJmEd97xhftEpO1SMERswAAfbvrb38K8eY0fd+aZ8O678Ic/+Kf1YnDYYdCpk5qTRNo6BUMeXHih79Nw3nkN//7hh+Gee+AnP/GhosVi663h0EO1qJ5IW6dgyIMePTwUJk2C55+v+7s1a/yO4nOfazw44pRMwpIl8PLLcVciInFRMOTJ974HvXptvlTGWWfBW295E1KHDvHV15ijj/YRVmpOEmm7FAx5stVWcPnlMH26jz4CeOIJ+NOf/E5h2LB462vMTjtBIuHNXSLSNikY8ujkk70P4fzzfRTSaafBHnt4H0QxSyZh1ixYtizuSkQkDgqGPGrf3rf/XLTIP4WvXOlNSJ06xV1Z0zKL6k2aFG8dIhIPBUOejR0Lhxziezj/8IcwfHjcFW3ZkCG+Q536GUTaJgVDnpnBrbf6hj5bmvRWTJJJ+Mc/4L334q5ERApNwVAAgwbBVVf5NqClIpn0GdxPPhl3JdIcmnciUVIwSIMSCd+6VKOTit/06T535oUX4q5EyoWCQRrUvr3PaXjySd9ASIpTba1v8vT22/DUU3FXI+VCwSCNSibh/fe9r0GK069/7XcKnTr5vt0iUVAwSKMOPdQn6ml0UnFatcrnxBx6qC/3/p//tJ2+hhdf9GXuV6+Ou5LyFEkwmNlYM3vVzBaa2WYrAJnZKDObZWYbzWx8vd+dYmY16a9ToqhHotGli6+4OmlS23nDKSXnnAMffQQ33eQ7AL75pq9z1RY88ADMmAF/+UvclZSnnIPBzNoDNwNHAEOBKjMbWu+wpcCpwN31XrsDcDEwEhgBXGxm2+dak0QnmYTly30mtBSP55+HO+/0cBg82AcLQNtpTkql/PHuu5s+TlonijuGEcDCEMLiEMJ64B4gmX1ACGFJCGE2sKneaw8HpoQQ1oQQ3gGmAGMjqEkictRR0K6dmpOKyYYN8O1vw2c/Cxdc4M/ttZff4bWFYKithWnTfJn4adN88qhEK4pg6AVkr6qzPP1cvl8rBbDjjr4ftIatFo9f/QrmzvXHrbby5yoqfFZ9WwiGOXPgww99DTLwvU0kWlEEgzXwXHNbpJv9WjM7zcyqzax6tXqcCiqZ9P0ZXnst7kpkxQqfQX/UUXDMMXV/l0j4CKVPPomltILJNCOdeCIccID6GfIhimBYDvTJ+rk38HrUrw0h3B5CqAwhVPbo0aNVhUrrZBbVU3NS/H7wA9i40e8W6kskvJmp3Ce6pVI++bJ/f6iq8juIOXPirqq8RBEMM4BBZtbPzDoCE4Dmrss5GTjMzLZPdzofln5OisiAAb5cuIIhXk8/DffdBz/+MfTrt/nvR470x3JvTkqlfBSWGRx/vE/G1F1DtHIOhhDCRuAM/A19HnBfCGGumV1mZscAmNlwM1sOHA/cZmZz069dA1yOh8sM4LL0c1Jkkkn45z99a1IpvHXr4IwzYOBAH4nUkF12gV13Le9gWL0aamo8GMDvHEaP9mDQkOroRDKPIYTwRAhhtxDCgBDCFennLgohTEp/PyOE0DuE0DWE0D2EsEfWa38fQhiY/vpDFPVI9MaN89Egjz8edyVt07XX+hviTTdB586NH5dIlHcwZP62TDAATJzo/V/TpsVTUznSzGdplv3280+kGp1UeEuWwBVXwHHHweGHN31sIgFLl/qmUOUolfIRWJWVnz537LG+JIiak6KjYJBmadfOR8FMnlz+o16Kzfe+5//8r7tuy8dmJrqV66fnVAr22efTYboA22wDX/qS97/U1sZXWzlRMEizJZOwdi0880zclbQdjz3mS5JcdBH06bPl44cNgw4dyrM5aeNGX2I8uxkpo6oK3nhDCz5GRcEgzfbFL0K3bhqdVCgffwxnngm77w7f/37zXtO5s4dDOQbD7Nm+NtTnP7/57770Jf9vU0tkREPBIM3WqZPvYf3oo7Cp/uImErmf/cw7VW+5BTp2bP7rEglfYG7jxvzVFofMxLaG7hi6dPEBEg884CO4JDcKBmmRceP8ln369LgrKW81NXD11XDSSfCFL7TstYmEf7Iut0lfqRTsvLMPyW3IxIm+R7k2LMqdgkFa5MgjfVSIRiflTwg+Z6FzZx+m2lLlutLq1KmfTmxryOjRvraXRiflTsEgLbLddnDwwepnyKcHHvBZzpdfDjvt1PLX9+3rE7/KKRhWrfJmtYb6FzI6dPCZ0JMm+SJ70noKBmmxZBLmz4cFC+KupPx8+KF3NO+zjy+t3Rpm5TfRran+hWxVVd5pP6m5i/JIgxQM0mKZVT111xC9yy7zFVRvucWb7ForkYBXXy2fJUxSKb8j2Hffpo874ADo3Vujk3KlYJAW++xn/ROtgiFac+f6JLavfnXLn4y3JNPPUC6DBKZO9VBoajkQ8ImAEyb4RMy33y5MbeVIwSCtMm6c/8/65ptxV1IeQoDvfMfH4l91Ve7nq6z0N8lyaE5avx6qq5vuX8g2caIP1X3ggfzWVc4UDNIqyaS/mT36aNyVlIe774bnnvO5C1FsN9KtG+y5Z3kEw0sv+TIszb2L2mcf3wdbo5NaT8EgrfK5z3mTkpqTcvfee3D22b4159e/Ht15EwlfM6nUJyNOneqPzQ0GM++Efu4576+RllMwSKuYeSf0lCm+fpK03k9+4k1yt9zim85EJZGAd98t/dFjqZR3KPfu3fzXVFX5He199+WvrnKmYJBWSyb9Fn/KlLgrKV0vvAA33wynn153KekolMtEt1Sq+f0LGbvt5p3VGp3UOgoGabVRo3zCm5qTWmfTJp+r0L07/PSn0Z9/8GDYdtvSDoYVK3x/idaM0qqq8k7rmpro6yp3CgZptQ4dfImMRx8tvwXbCuEPf/A37Wuuge23j/787dr5PtClHAzNndjWkAkTvMnznnuiraktUDBITsaN8/HimQ5CaZ6334Zzz4UDD4STT87fdRIJePnl0l0iIpXyVX2HDWv5a3v3hoMO0n7QraFgkJyMHetLQqs5qWV+/GPvGL7llsYXhYtCIuFNVtXV+btGPqVS3vfSkmXHs1VVwbx5vpeDNJ+CQXLSrRsccogHQzF8KtuwwWe9FvPEu+nT4Te/8U149torv9caMcIfS7E5ad06mDkzt1ng48f70iLqhG4ZBYPkLJmERYvglVfiq6G2Fu68E4YM8buYAQN83aFia0KprfURSDvtBJdckv/rde/uI3RKMRhmzfJZz7kEw447wpgx3s9Q6vM5CknBIDmLc1G9EHzpg7339rb6bt08IA47DC6+GAYNgttvL57O8dtu8ze8X/7SN7EvhEwHdDHc0bVELh3P2SZO9JFNmfPJlikYJGe77OKzdgsZDCHAk096+/P48f5p8L77/E33K1/xsPj3v/3O4Zvf9CabuJu7Vq3yvoVDDoETTyzcdRMJv/Z//1u4a0YhlfK9JXbeObfzJJO++J6WyGg+BYNEYtw4bzt//fX8X+u553y0yZFH+rLSf/yjj7w5/ngfopnx+c/DP//pu82F4DWOGhVfs8qPfuRbbt58c347nOsrxYluIXy6Y1uuunWDo4+Gv/61eO4ci52CQSKRTPpjPjdImT7dm4i+8AXfzeuWW3zPgVNOaXzvAjOvbc4cuPVWWLjQ32zGjy/sUhH//Cf86U++JtKQIYW7LvjdUpcupRUMy5b5h4woggF8dNKbb8Kzz0ZzvnKnYJBIDB3qzTb5aE6aPds/7Y8c6UtIXHutv8GffnrzhzFWVHiTUk0NXHqpj1waOtSXul61Kvqas23Y4DOcd90VLrwwv9dqSIcO3uRWSsGQ6Q9o6VIYjTniCO/T0eik5lEwSCQyn8yffRY++CCacy5Y4J/09tkH/v53H2W0eLF/6u7SpXXn3HpruOgiD5ZvftM7pgcOzO8Iphtv9DuWG26Arl3zc40tSSQ8VNeti+f6LZVK+b/jvfeO5nydO8OXvwwPPeTre0nTIgkGMxtrZq+a2UIzO6+B33cys3vTv59mZn3Tz/c1s4/N7MX0161R1CPxSCZ9eOFTT+V2nqVLffnpoUO9aercc73p6Cc/8fbiKHzmM97WP3cuHH74pyOYbrst2nboFSv83Ece+WlzWxwSCf9388IL8dXQEqmUD2jo0CG6c06cCO+/D088Ed05y1YIIacvoD2wCOgPdAReAobWO+bbwK3p7ycA96a/7wvMaek199tvvyDFZ8OGELp3D+Gkk1r3+pUrQ/jud0Po2NG/zjzTnyuEqVNDOPDAECCEIUNCeOihEDZtyv28J54YQqdOISxcmPu5crFihf9t110Xbx3N8dFHIVRUhHDuudGed8OGEHr2DGH8+GjPW0qA6tCM99go7hhGAAtDCItDCOuBe4D6n42SwB3p7+8HRpsVclyGFEJFhY/+ePxxb1dvrrffhvPOg/79vUP55JO9L+CGG3wiWCHsvz88/7yPYAI49lgf+ZTL2PdnnoF774Xzz/f+lzjtsgv06VMa/QwzZ/pdW1T9CxkVFXDCCfDYY37nII2LIhh6Acuyfl6efq7BY0IIG4H3gO7p3/UzsxfM7DkzO6ixi5jZaWZWbWbVq1evjqBsyYdk0tcAev75LR/7/vveEdy/P/z85/5mPG+eLxex6675r7W+TD/Jyy97k9KiRf7mdNxxLR/BtG6dd2wPGOBNYcUgkSiNYMiEcWaYbZSqqryPQWt7NS2KYGjok3/9aUSNHbMS2DWEMAz4AXC3mTU4HzSEcHsIoTKEUNkjik1xJS/GjPGOvqb+x/voI19qun9/XxZi9GgfeXTXXd7OH7eKCjjtNO+gvuwyePpp7+/49rebP4LpF7/wobQ33uj/PIpBIuGT3FaujLuSpqVSHqg9e0Z/7v339y1pNTqpaVEEw3KgT9bPvYH605z+d4yZVQDbAmtCCOtCCG8DhBBm4n0Vu0VQk8Ska1cPh4ZmGa9f7x2+Awf6ZK/KSpgxAx580DeuLzZdu3qH96JF8K1v+Z3MgAF+l9PUCKYlS3zjnWOP9WGSxSLzCXzatHjraEqUE9saYub7NEyZAmp4aFwUwTADGGRm/cysI965XH+a0yTglPT344FnQwjBzHqYWXsAM+sPDAIWR1CTxCiZ9JFFL73kP2/c6JvS7LYbnHGGB8Nzz/nopai3s8yHnj3hppt8kcAjjvC7nIEDfcJcQ30p3/++vwFdf33BS23SsGE+yqeYm5OWLPG7sqj7F7JNnOiLGd5/f/6uUepyDoZ0n8EZwGRgHnBfCGGumV1mZunl1fgd0N3MFuJNRpkhraOA2Wb2Et4p/a0Qwppca5J4HXWUvzE+/LB3vu65J3z1q77S5VNPeSiMGhV3lS03aJAvq5BKecidfrrPKs4suQHe8f7II36nEUc/SVO6dPE5IcUcDFEtnNeUvfbypkGtndSE5gxdKrYvDVctfgccEIKZD5HcY48QHnwwmuGfxWLTphAeeSSE3Xf3v/Hznw/hmWdC6NfPh7uuWxd3hQ377ndD2GorH7pZjM44I4SuXfNf3+WX+7+3pUvze51iQwGHq4ps5jvfgX33hT//2ZuUjj22sAvH5ZuZLzc+e7bPnl682DvRX3vN+1Fau+NYviUS3vk/Z07clTRs6lRf+qSxta+iMmGCP2o/6IYpGCQvqqp8O8mTToL27eOuJn8qKuAb3/ARTFdeCVdd5ctqF6tiXml17Vr/EJHPZqSMgQN9ZrWakxqmYBCJQNeuPpGtWOYsNKZfP+jRoziDobraO4ULEQzgH15eeMGHFUtdCgaRNsSseCe6TZ3qj/mY2NaQE0/0fx66a9icgkGkjUkk/FPyO+/EXUldqRQMHuz7VBfCLrv43h5/+UvpbXuabwoGkTYm84l8+vR468gWggdDoZqRMqqqfLmTUll1tlAUDCJtzPDh3oRSTM1JixbBW28VPhiOO84n/WmJjLoUDCJtTLduPumwmIIh079Q6GDYYQffj+Pee2HTpsJeu5gpGETaoETC10wqljfDVMq33hw6tPDXrqqC5cvhX/8q/LWLlYJBpA1KJLzzuaYm7kpcKuUT2+KY83LMMbDVVhqdlE3BINIGFdNEtw8+8D0wCt2MlLH11h4Of/1ryzaYKmcKBpE2aMgQb7ophmCYPt2btOIKBvDmpLffhr/9Lb4aiomCQaQNatfOm26KIRjyuWNbcx1+OGy3nUYnZSgYRNqoRMIXAVy7Nt46UinvdN5uu/hq6NTJh64+/LAvMtjWKRhE2qhEwptwqqvjq2HTJr9ribMZKWPiRN+Z7/HH464kfgoGkTZq5Eh/jLM5acECWLOmOILh4INh5501OgkUDCJtVvfuvitdnMGQ6V/I51aezdW+PZxwAjzxBLz7btzVxEvBINKGZVZajWsRuVTK+xYGD47n+vVVVcG6dfDQQ3FXEi8Fg0gblkjAG2/A0qXxXD+V8hraFck70YgR0L+/mpOK5F+HiMQhzolu770Hc+cWRzNShplv+/nMM7BqVdzVbK62tjDXUTCItGF77QVdusQTDNOmeRNWMXQ8Z5s40UdL/fWvcVfyqU2b4IYboLISPv44/9dTMIi0YR06+JtNHMGQSvkn9BEjCn/tpuyxhwdmsTQnLVzoI6a+/33o1asw804UDCJtXCIBs2Z5p2shpVK+/Pc22xT2us1RVeVLgS9ZEl8NmzbB9dfD3nvDnDlwxx3w6KOw4475v7aCQaSNSyRg/Xp48cXCXTMzsa2Y+heyTZjgj/fcE8/1a2r8LuGss2D0aO+LOflkv8MqBAWDSBsXRwf0vHne+Vxs/QsZ/fr5P5dCNyfV1vpdwuc+53cJf/oTTJrk+1MXkoJBpI3bZRfo06ewwZCZ2FaswQDenDR7NrzySmGu19Bdwv/9X+HuErIpGETkfxPdCiWV+nTmdbE64QSfX5Hvu4baWrjuOu9LmDs3vruEbAoGESGR8I7WN94ozPWmTvW7hTg+DTfXTjvBIYd4MORrZviCBX6X8IMfwJgxfncS111CtkiCwczGmtmrZrbQzM5r4PedzOze9O+nmVnfrN+dn37+VTM7PIp6RKRlMv0M06bl/1pr1sD8+cXdjJRRVQWLFkW/Am1tLfzyl96X8MorcOed8MgjvohfMcg5GMysPXAzcAQwFKgys/pben8NeCeEMBC4Drg6/dqhwARgD2AscEv6fCJSQMOGQUVFYZqTMuFTCsHw5S9Dx47RbuCzYAGMGgVnn+13CXPnwle+Ev9dQrYo7hhGAAtDCItDCOuBe4BkvWOSwB3p7+8HRpuZpZ+/J4SwLoTwGrAwfT4RKaAuXWCffQoTDFOnetv98OH5v1auttsOjjgC7r039+Uosu8S5s2DP/+5uO4SskURDL2AZVk/L08/1+AxIYSNwHtA92a+FgAzO83Mqs2sevXq1RGULSLZEgmYMQM2bszvdVIpf3Pceuv8XicqVVWwciU8/3zrz/Hqq5/eJRx2mN8lnHRScd0lZIsiGBr60+p31TR2THNe60+GcHsIoTKEUNmjR48WligiW5JI+HILc+dUM5jxAAAKxUlEQVTm7xq1td6UVArNSBlHHw1du7ZudFJtLfziF343lrlLePjh4rxLyBZFMCwH+mT93Bt4vbFjzKwC2BZY08zXikgBFGKi29y5vn1mKQXDVlvBuHFw//0+Q7y5Xn0VDjoIfvjD0rhLyBZFMMwABplZPzPriHcmT6p3zCTglPT344FnQwgh/fyE9KilfsAgYHoENYlIC/Xv7+vw5DMYpk71x1IKBvDmpHfegcmTt3xsbS1ce63fJcyfXzp3CdlyDoZ0n8EZwGRgHnBfCGGumV1mZsekD/sd0N3MFgI/AM5Lv3YucB/wCvAU8J0QQoFWHBeRbGb5n+iWSkHPnh5CpWTMGNhhhy03J82fDwceCOecA4cf7kNRS+UuIVtFFCcJITwBPFHvuYuyvv8EOL6R114BXBFFHSKSm0QCHnvMPx1vv33050+lin9iW0M6doTx4/3T/9q13ueQLTN7+cIL/Xd33eV3GaX2d2Zo5rOI/E+mn2F6Hhp033rL1wMqtWakjIkT4aOPfOnrbNl3CWPHel/CxImlGwqgYBCRLMOH+xtaPpqTMgvnFetS21ty0EG+UU6mOam2Fq65xvsSFizwSXAPPeRLaZQ6BYOI/M822/gOZvkKhooK3zGuFLVrByeeCE8+6X/LAQfAj37kE+Dmzi3tpqP6FAwiUkci4XMNNm2K9ryplH+67tIl2vMWUlUVbNjgdz01NX6X8OCD5XGXkE3BICJ1JBLe+VxTE905N270fotS7V/I2G8/H210/PE+4qic7hKyRTIqSUTKR/ZEt8GDoznn7NnecVuq/QsZZvDUU3FXkX+6YxCROnbf3fsaouxnKIUd2+RTCgYRqaNdOxgxIvpg2Hln2HXX6M4p+aNgEJHNJBLe/LN2bTTnK9WJbW2VgkFENpNI+KikKHYuW7UKFi8u/f6FtkTBICKbGTnSH6PY6lP9C6VHwSAim9lxRxg4MJp+hlQKOnSAfffN/VxSGAoGEWlQIuFv6qHBrbOaL5XyUOjcOZq6JP8UDCLSoEQC3ngDli3b8rGNWb/etwtV/0JpUTCISIOi2NHtpZfgk0/Uv1BqFAwi0qC99/bmn1yCQR3PpUnBICIN6tDBV0LNNRj69IHevaOrS/JPwSAijUokYNYsWLeuda+fOlV3C6VIwSAijUokPBReeqnlr339dVi6VMFQihQMItKoXDqg1b9QuhQMItKoXr28f6C1wdCpEwwbFn1dkl8KBhFpUiLRumCYOtU7rzt2jL4myS8Fg4g0KZGA117zxfCaa906mDlTzUilSsEgIk3K9DO0ZEG9F17wWc8KhtKkYBCRJu27L1RUtKw5aepUf1QwlCYFg4g0qUsX2GeflgVDKgV9+/qubVJ6FAwiskWJBEyfDrW1zTs+s2OblKacgsHMdjCzKWZWk37cvpHjTkkfU2Nmp2Q9/w8ze9XMXkx/9cylHhHJj0TCt/mcO3fLxy5bBitWKBhKWa53DOcBz4QQBgHPpH+uw8x2AC4GRgIjgIvrBchJIYR90l9v5liPiORBSya6ZfoXtNR26co1GJLAHenv7wDGNXDM4cCUEMKaEMI7wBRgbI7XFZEC6t/fd3VrTjCkUt4vsffe+a9L8iPXYPhMCGElQPqxoaagXkD2Vh/L089l/CHdjPQTM7Mc6xGRPDBr/kS3VAqGD/fVWaU0bTEYzOxvZjanga9kM6/R0Jt9ZrPAk0IIewEHpb/+r4k6TjOzajOrXr16dTMvLSJRSSRg3jx4993Gj/n4Y5/DoP6F0rbFYAghHBpC2LOBr0eAVWa2M0D6saE+guVAn6yfewOvp8+9Iv34AXA33gfRWB23hxAqQwiVPXr0aO7fJyIRyfQzTJ/e+DEzZ8KGDepfKHW5NiVNAjKjjE4BHmngmMnAYWa2fbrT+TBgsplVmNmOAGbWATgKmJNjPSKSJ8OHe5NSU81JmRVVMyEipSnXYLgKGGNmNcCY9M+YWaWZ/RYghLAGuByYkf66LP1cJzwgZgMvAiuA3+RYj4jkyTbbwB57bDkYBgyAnhp4XtIqcnlxCOFtYHQDz1cDX8/6+ffA7+sdsxbYL5fri0hhJRLwwAMQgt89ZAvBg2HMmHhqk+ho5rOINNvIkfDOO1BTs/nvliyBN95Qx3M5UDCISLM1NdFNO7aVDwWDiDTb7rtDt26NB0PXrrDnnoWvS6KlYBCRZmvfHkaMaDwYRo70JbqltCkYRKRFEgmYPdsX1ctYuxZefFHNSOVCwSAiLZJI+PLbM2d++lx1tT+nYCgPCgYRaZGRI/0xuzlJE9vKi4JBRFqkRw+fxFY/GAYPhu7d46tLoqNgEJEWSyQ8DELwr6lT1YxUThQMItJiiYRPZlu2DBYtgrfeUjCUEw0sE5EWy57otm6df69gKB8KBhFpsb33hs6dPRg++cQX2Bs6NO6qJCoKBhFpsY4dYb/9PBg++shHKrVvH3dVEhX1MYhIqyQSPpfh5ZfVjFRuFAwi0iqJBKxfD5s2KRjKjYJBRFolezKbJraVF/UxiEir9O4NvXrBttvCdtvFXY1EScEgIq12zTU+OknKi4JBRFqtqiruCiQf1McgIiJ1KBhERKQOBYOIiNShYBARkToUDCIiUoeCQURE6lAwiIhIHQoGERGpw0IIcdfQYma2Gvhv3HW00I7AW3EXUWD6m9sG/c2l47MhhB5bOqgkg6EUmVl1CKEy7joKSX9z26C/ufyoKUlEROpQMIiISB0KhsK5Pe4CYqC/uW3Q31xm1McgIiJ16I5BRETqUDAUgJltZ2b3m9l8M5tnZmW/Q66ZnWVmc81sjpn9xczKbjsXM/u9mb1pZnOyntvBzKaYWU36cfs4a4xaI3/zNen/tmeb2UNmVlb7uTX0N2f97odmFsxsxzhqyxcFQ2HcADwVQhgCfA6YF3M9eWVmvYAzgcoQwp5Ae2BCvFXlxR+BsfWeOw94JoQwCHgm/XM5+SOb/81TgD1DCHsDC4DzC11Unv2Rzf9mzKwPMAZYWuiC8k3BkGdmtg0wCvgdQAhhfQjh3XirKogKoIuZVQBbAa/HXE/kQgjPA2vqPZ0E7kh/fwcwrqBF5VlDf3MI4ekQwsb0j/8Behe8sDxq5N8zwHXAj4Cy66hVMORff2A18Acze8HMfmtmXeMuKp9CCCuAa/FPUiuB90IIT8dbVcF8JoSwEiD92DPmegrtq8CTcReRb2Z2DLAihPBS3LXkg4Ih/yqAfYFfhxCGAWspv+aFOtLt6kmgH7AL0NXMvhJvVZJvZnYBsBG4K+5a8snMtgIuAC6Ku5Z8UTDk33JgeQhhWvrn+/GgKGeHAq+FEFaHEDYADwKfj7mmQlllZjsDpB/fjLmegjCzU4CjgJNC+Y+BH4B/6HnJzJbgTWezzGynWKuKkIIhz0IIbwDLzGxw+qnRwCsxllQIS4GEmW1lZob/zWXd4Z5lEnBK+vtTgEdirKUgzGwscC5wTAjho7jrybcQwsshhJ4hhL4hhL74h7990/+vlwUFQ2F8F7jLzGYD+wBXxlxPXqXvju4HZgEv4/+dld1MUTP7C5ACBpvZcjP7GnAVMMbMavARK1fFWWPUGvmbbwK6AVPM7EUzuzXWIiPWyN9c1jTzWURE6tAdg4iI1KFgEBGROhQMIiJSh4JBRETqUDCIiEgdCgYREalDwSAiInUoGEREpI7/DzUkRK+WXiV4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x330096def0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 绘制不同聚类数目的模型性能， 找到最佳模型\n",
    "plt.plot(Ks, np.array(CH_scores), 'b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最佳值还是在边界，在5左侧再试一次。。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K-means begin with 2 clusters\n",
      "CH_score: 0.6375833065268858\n",
      "K-means begin with 3 clusters\n",
      "CH_score: -0.028437240892212824\n",
      "K-means begin with 4 clusters\n",
      "CH_score: 0.2654092150757131\n",
      "K-means begin with 5 clusters\n",
      "CH_score: 0.04129695217077523\n",
      "K-means begin with 6 clusters\n",
      "CH_score: 0.11633798272466671\n",
      "K-means begin with 7 clusters\n",
      "CH_score: 0.07346162486032493\n",
      "K-means begin with 8 clusters\n",
      "CH_score: 0.2023301543560187\n",
      "K-means begin with 9 clusters\n",
      "CH_score: 0.11293376045929865\n",
      "K-means begin with 10 clusters\n",
      "CH_score: 0.04793811774937678\n"
     ]
    }
   ],
   "source": [
    "# 调整超参数（聚类数目K）的搜索范围\n",
    "CH_scores = []\n",
    "Ks = [2,3,4,5,6,7,8,9,10]\n",
    "for K in Ks:\n",
    "    ch = K_cluster_analysis(K, eventContMatrix.toarray())\n",
    "    CH_scores.append(ch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x3300a04470>]"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XmYVPW17vHvYlJBHEGjgIKKAwJSTV+ikpgY9IhDQK/GA0k8mmgccc450ZMTB4zzxCBRUWOch5hBVKIxJnqiN0ZbAyogiqKAaMRZcYCm1/1jVcWy7aaL7qraVbvez/P0013Vm70X3dVv7f3bv8HcHRERSZdOSRcgIiLFp3AXEUkhhbuISAop3EVEUkjhLiKSQgp3EZEUUriLiKSQwl1EJIUU7iIiKdQlqQP36tXL+/fvn9ThRUSq0lNPPfWWu/dua7vEwr1///40NDQkdXgRkapkZq8Wsp2aZUREUkjhLiKSQgp3EZEUUriLiKSQwl1EJIUU7iIiKaRwFxFJoaoL98ceg9NPB60OKCLSuqoL96efhgsugKVLk65ERKRyVV24ZzLx+R//SLYOEZFKVnXhvtNOYKZwFxFZnaoL9549YeBAhbuIyOpUXbhDNM08/XTSVYiIVK6qDfdXX4V33km6EhGRylSV4V5XF59nzUq2DhGRSlWV4Z7rMaOmGRGRlhUU7mY22szmm9kCMzutlW0ONrO5ZjbHzG4tbplf1KsX9O2rm6oiIq1pcyUmM+sMTAP2BJYAT5rZDHefm7fNQOB0YKS7v2tmm5Sq4Jy6OoW7iEhrCjlzHwEscPeX3X0FcDswttk2PwKmufu7AO7+ZnHL/LJMBubPh+XLS30kEZHqU0i49wEW5z1ekn0u37bAtmb2mJk9bmajW9qRmR1pZg1m1rBs2bL2VZyVyUBTEzzzTId2IyKSSoWEu7XwXPNpu7oAA4FvAuOBa81sgy/9I/fp7l7v7vW9e7e5ePdqaRoCEZHWFRLuS4B+eY/7As2n7VoC3O3uK919ITCfCPuS6dcPNt5Y4S4i0pJCwv1JYKCZDTCzbsA4YEazbX4P7A5gZr2IZpqXi1loc2Zx9q5wFxH5sjbD3d0bgQnAA8A84E53n2NmE81sTHazB4C3zWwu8BfgP9397VIVnZPJwLPPwsqVpT6SiEh1abMrJIC7zwRmNnvujLyvHTgl+1E2mQysWAFz58ZskSIiEqpyhGpObhoCNc2IiHxRVYf7wIHQo4emIRARaa6qw71Tp2iO0Zm7iMgXVXW4QzTNzJoVA5pERCRUfbhnMvDRR/DSS0lXIiJSOVIR7qB2dxGRfFUf7jvuCF27qt1dRCRf1Yd7t24weLDCXUQkX9WHO3w+DYE3n85MRKRGpSbcly2D115LuhIRkcqQmnAHNc2IiOSkItx32ilmiVS4i4iEVIT7uuvCttuqO6SISE4qwh00t7uISL7UhHtdHSxaBG+XfBZ5EZHKl5pwz91UnTUr2TpERCpB6sJd7e4iIikK9403jkWz1e4uIpKicIdod1e4i4ikLNwzGZg/H5YvT7oSEZFkFRTuZjbazOab2QIzO62F7x9mZsvMbFb244jil9q2TCbml5k9O4mji4hUjjbD3cw6A9OAvYFBwHgzG9TCpne4+7Dsx7VFrrMgmoZARCQUcuY+Aljg7i+7+wrgdmBsactqn759oVcvhbuISCHh3gdYnPd4Sfa55g40s2fM7C4z61eU6taQWZy9qzukiNS6QsLdWniu+czp9wD93X0o8CfghhZ3ZHakmTWYWcOyZcvWrNICZTLw3HOwYkVJdi8iUhUKCfclQP6ZeF9gaf4G7v62u3+WfXgNMLylHbn7dHevd/f63r17t6feNtXVwcqVMHduSXYvIlIVCgn3J4GBZjbAzLoB44AZ+RuY2WZ5D8cA84pX4prRTVURkQLC3d0bgQnAA0Ro3+nuc8xsopmNyW52gpnNMbPZwAnAYaUquC3bbBNTAKvdXURqmXlCC4/W19d7Q0NDSfb9ta/F50cfLcnuRUQSY2ZPuXt9W9ulaoRqTl1dDGRqakq6EhGRZKQy3DMZ+OgjWLAg6UpERJKR2nAHtbuLSO1KZbgPGgRdu6rHjIjUrlSGe7duMGSIwl1Ealcqwx0+n4Ygoc5AIiKJSnW4v/02LFmSdCUiIuWX6nAHNc2ISG1KbbjvtFPMEqlwF5FalNpw79EDtttO3SFFpDalNtwhmmZ05i4itSjV4V5XB4sXw1tvJV2JiEh5pTrcdVNVRGqVwl1EJIVSHe4bbQRbbKFwF5Hak+pwh2h3V7iLSK1JfbhnMvDCCzEFsIhIraiJcHePxTtERGpFTYQ7qGlGRGpL6sO9Tx/o3VvhLiK1JfXhbvb59L8iIrWioHA3s9FmNt/MFpjZaavZ7iAzczNrc2XucspkYM4cWLEi6UpERMqjzXA3s87ANGBvYBAw3swGtbBdT+AE4O/FLrKj6upg5coIeBGRWlDImfsIYIG7v+zuK4DbgbEtbHcOcBHwaRHrKwotmC0itaaQcO8DLM57vCT73L+YWQbo5+73FrG2otl6a+jZUzdVRaR2FBLu1sJz/1qZ1Mw6AZcDp7a5I7MjzazBzBqWLVtWeJUd1KlTLN6hcBeRWlFIuC8B+uU97gsszXvcExgMPGxmrwA7AzNauqnq7tPdvd7d63v37t3+qtuhri4GMq1aVdbDiogkopBwfxIYaGYDzKwbMA6Ykfumu7/v7r3cvb+79wceB8a4e0NJKm6nTAaWL4cXX0y6EhGR0msz3N29EZgAPADMA+509zlmNtHMxpS6wGLRSFURqSVdCtnI3WcCM5s9d0Yr236z42UV36BB0K1bhPv48UlXIyJSWqkfoZrTtSsMGaIzdxGpDTUT7vD5NATubW8rIlLNai7c33knFs0WEUmzmgr3urr4rKYZEUm7mgr3oUNjQJOmIRCRtKupcO/eHbbbTmfuIpJ+NRXuEO3uCncRSbuaC/e6OliyBMo4tY2ISNnVXLhrpKqI1IKaC/dhw+Kzwl1E0qzmwn2jjWDLLRXuIpJuNRfuEO3uCncRSbOaDPdMBl54AT78MOlKRERKo2bDHWLxDhGRNKrJcNc0BCKSdjUZ7pttBptsomkIRCS9ajLczTRSVUTSrSbDHSLc58yBzz5LuhIRkeKr2XCvq4PGxgh4EZG0qdlwz/WYUbu7iKRRzYb7VltBz55qdxeRdCoo3M1stJnNN7MFZnZaC98/2syeNbNZZvaomQ0qfqnF1alTzDOjcBeRNGoz3M2sMzAN2BsYBIxvIbxvdfch7j4MuAi4rOiVlkBdXQxkWrUq6UpERIqrkDP3EcACd3/Z3VcAtwNj8zdw9w/yHvYAvHgllk4mAx9/HFMRiIikSSHh3gdYnPd4Sfa5LzCz48zsJeLM/YSWdmRmR5pZg5k1LKuA1TI0t7uIpFUh4W4tPPelM3N3n+buWwM/Af6npR25+3R3r3f3+t69e69ZpSWwww6w1loKdxFJn0LCfQnQL+9xX2Dpara/Hdi/I0WVS9euMGSIukOKSPoUEu5PAgPNbICZdQPGATPyNzCzgXkP9wVeLF6JpZWbhsCr4i6BiEhh2gx3d28EJgAPAPOAO919jplNNLMx2c0mmNkcM5sFnAIcWrKKiyyTgXffhUWLkq5ERKR4uhSykbvPBGY2e+6MvK9PLHJdZZM//e+WWyZbi4hIsdTsCNWcIUNiQJPa3UUkTWo+3Lt3h+23V48ZEUmXmg930NzuIpI+Cnei3f211+DNN5OupDTUE0ik9ijcSfdI1bvvhg03hPnzk65ERMpJ4U7MDgnpC3d3OPdceP99OPXUpKsRkXJSuBNntgMGpC/cH38cnnwymp3uuw8eeCDpikSkXBTuWZlM+rpDTpoE668PDz4IW28Np5wSSwuKSPop3LMyGViwAD74oO1tq8HixfCb38CPfgQbbQSXXAJz58LVVyddmYiUg8I9K3dTdfbsZOsolmnTos19woR4PHYs7L47nHFGTLcgIummcM/Kn4ag2i1fDtOnwwEHfD6lglk007z3Hpx9drL1iUjpKdyzNtsMNt00He3uN98cZ+cnnfTF54cOhSOOiLP6559PpjYRKQ+Fe540jFRtaooz9OHDYeTIL3//nHNiygV1jRRJN4V7nkwmbjp+9lnSlbTfgw/GWfmJJ0ZTTHObbAI/+xnMnAn331/++kSkPBTueerqoqvgc88lXUn7TZ4MX/kKHHxw69uccAJss426RoqkmcI9T67HTLW2uz//PPzhD3DssbE2bGu6dYuukfPmwVVXla8+ESkfhXueAQNgvfWqt919ypQI7qOOanvbMWPgW9+CM8+Ed94pfW0iUl4K9zydOlXvTdV334UbboDvfS/a1dtiBpdfrq6RImmlcG8mk4mBTKtWJV3Jmrn2Wvj447iRWqihQ2ME67Rp0UQjIumhcG8mk4FPPqmuKXIbG2HqVPjmN2Gnndbs355zDvTooa6RImlTULib2Wgzm29mC8zstBa+f4qZzTWzZ8zsITOr2qWmq3Fu99//PuaSaT5oqRC9e8eUBH/4Q3yISDq0Ge5m1hmYBuwNDALGm9mgZpv9A6h396HAXcBFxS60XHbYAdZeu7rCfdIk2Gor2G+/9v3744+PrpGnngorVxa3NhFJRiFn7iOABe7+sruvAG4HxuZv4O5/cfePsw8fB/oWt8zy6dIFhgypnu6QDQ3w2GMR0J07t28f3brBpZeqa6SsOfd4/VXzwL+0KiTc+wCL8x4vyT7XmsOBqr7Az/WYqYa1RydPhp494Yc/7Nh+vv1tGDVKXSOlcO5xA/9rX4MRI+DZZ5OuSPIVEu4tDGKnxdgzs+8D9cDFrXz/SDNrMLOGZcuWFV5lmWUy0UXw1VeTrmT1li6FO+6IYF9vvY7tK9c18v334ayzilKepJh73OOZOhX+/d/hjTegvh4uuyzmN5LkFRLuS4B+eY/7Akubb2RmewA/Bca4e4sXae4+3d3r3b2+d+/e7am3LKpl+t8rr4yeMscfX5z9DRkCRx4Jv/iFukZK69zh5JNj0NxJJ8Ftt8VZ++jRcd9mjz1g0aKkq5RCwv1JYKCZDTCzbsA4YEb+BmaWAa4mgv3N4pdZXkOGRPt1Jbe7f/pptI9/+9uxhF6xTJwI666rrpHSMvd4bUyeHE0yl10WV32bbBK9tq65Bp54IsZQ3HJLdTRtplWb4e7ujcAE4AFgHnCnu88xs4lmNia72cXAusCvzWyWmc1oZXdVYZ11YPvtK/vM/dZb4a232tf9cXXUNVJa4w4//nE0351wQnzOn3nULNYLmD0bBg2C738fxo/XPZykmCf01lpfX+8NDQ2JHLsQhxwCf/4zvPZa0pV8mXsMVjKDWbNantq3I1asgMGD4+rlmWega9fi7l+qjzv8539Gr6oJE6JJZnWvu8ZGuPDCuH+z6abwq19Fc410nJk95e71bW2nEaqtqKuLG5b//GfSlXzZww9HG2drc7Z3VK5r5PPPR7u+1DZ3+MlP4jVx3HFtBztEl+Kf/hT+9rdo5ttzz7jK/OST8tQsCvdWVfJI1UmToFcv+O53S3eM/faLM62zzoK33y7dcaSyucNpp8HFF8dU0lOnrtkJRX193Ls6/vhop6+vr8y/qTRSuLdi2LD4XGkvxJdegnvugaOPjpG0paKukeIOp58OF10ExxwDV1zRvivF7t3jbP/++2P20q9+Fc4/v/om56s2CvdWbLBBDOmvtHCfOjUueY85pvTHGjw45oa/8spYflBqh3s0q1x4YbwG2hvs+fbaK5oT998f/vu/4RvfgIULi1OvfJnCfTUymcrqDvnBB/DLX8agkc03L88x1TWy9rjD//xPnF3nxj10KlJSbLxxDLy78cYI+qFD4frr1WWyFBTuq5HJRDPI++8nXUn45S/hww/XbM72jurVK6YkuP/+WFRb0s09FlA/77yY6//KK4sX7Dlm0RvtmWdg+PAYYX3ggVDBg9arksJ9NXI3VWfPTrYOiPbJKVNg5Mi4KVVOxx0H224bC2pr1sj0co838nPPjf7qV11V/GDPt+WW8NBDcbP2vvti8KBOIIpH4b4alTQNwb33RvtksQctFSLXNXL+/LhEl3Q666xYvOXww+Hqq0sb7DmdO8fAqCefjAF0++4b95OWLy/9sdNO4b4aX/lKfFRCu/ukSbDFFnEzKgn77ht9ldU1Mp3OPjvur/zgBzB9enmCPd/QoRHwp54abyyZTExjIO2ncG9DJSyYPWtWDFyaMCF6yiQh1zXygw/i0l3SY+LEeNM+7LBYi7fcwZ6z9tpwySXRVPPpp7DrrvGm09iYTD3VTuHehkwmugF++mlyNUyZEn2FjzgiuRoAdtwx+tdfdRXMmZNsLVIc55wTb9aHHppssOfbffe42Tp+fLzpjBwJL76YdFXVpwJ+lZWtri5uZia1EMGbb8bseoceChtumEwN+c4+OxYHOeUUdV+rdueeG5PEHXIIXHdd+1fyKoUNNoCbbopuky++GIMKr75ar7k1oXBvQ9LTEFx1VUzkdcIJyRy/uVzXyD/+UT0bqtl550Vf9kMOiX7mlRTs+Q4+OE6sRo6Mq8b99ouFQaRtCvc2DBgA66+fTLh/9ln0M95775iCuFIce6y6Rlaz88+P0aff+15lB3tOnz4xzmLKlJipdciQmDteVk/h3gaz5G6q3nlnnKWUc9BSIbp1i0UaXngBpk1LuhpZExdeGEP/v/tduOGGyg/2nE6dYvKxp56Cfv3ggANi8NOHHyZdWeVSuBcgk4mBTOW8a+8e3R932AH+7d/Kd9xC7bNP1HX22bFoiFS+iy6KGR7Hj6+uYM83aBA8/ni8Qd1wQ6xr8OijSVdVmRTuBchkorfM/PnlO+Zjj0X/+lLN2d5RZnH2/uGHmjWyGlx8cczJPm5czOuSVJfaYujWLW4G/+//xuNvfCPCfsWKZOuqNAr3AiRxU3XSpOgdc8gh5TvmmlLXyOpw6aXwX/8VE87ddFN1B3u+kSPjivqww+I+ws47a/bSfAr3Amy/fQywKFe4v/IK/O53MSNf9+7lOWZ7nXVWdI08+WR1U6tEl10Ww/sPPhhuvjk9wZ7Ts2d04/zd72Dx4piIbMoUaGpKurLkKdwL0KVLDI8u1zQE06ZFs8dxx5XneB3Rq1cE/IMPxuRP1cw9Fh7fa6/o7lntVyOXXx7D+b/znRgrkbZgz7f//tFlctSoaMocPboy1z8uK3dP5GP48OFeTY46yn2DDdybmkp7nA8/dF9/ffeDDy7tcYppxQr37bZz33Zb988+S7qa9nn6afeRI93BvV8/d7P4eocd3M880/2555KucM1cfnnUf+CB8fupFU1N7ldf7d69u/uGG7pPmuT+8cdJV1VcQIMXkLEFnbmb2Wgzm29mC8zstBa+v5uZPW1mjWZ2UNHfgSpAJgPvvRdNJqV0ww0xf3wSsz+2V9eu1ds1ctmyWGlo+PCo/7rr4ne8dGn8XzbZJOZeGTw47jGcfXblt+tOnhzNZAceCLfdFr+fWmEWzZmzZsXf7EknwTbbxEpSn32WdHVl1lb6A52Bl4CtgG7AbGBQs236A0OBG4GDCnlXqbYz9yeeiDOhu+4q3TFWrYqz3xEjSn+FUGxNTe577RVXHcuWJV1N21audJ8yJa7GunRxP/lk93ffbXnbpUvdr7jCfbfdPj+j33FH97PPdp87t7x1t2XKlKjvgANq64y9NX/5S/zewL1vX/crr6zeq8scCjxzLyTcdwEeyHt8OnB6K9v+Kq3h/skn7p07u//0p6U7xn33xW/klltKd4xSmjMnfkbHHpt0Jav30EPugwfHz3rPPaPuQi1d6j516heDfvDgCPp580pXcyGmTo169t+/+gOsmJqa3P/0J/ddd42fz5Zbul9zTfW++RUz3A8Crs17fAhwRSvbpjbc3eOPeJ99Srf/Pfd033zz6v7DnDDBvVMn92efTbqSL1u4MNqgwX3AAPff/75jV0ivvRZnyl//+udBP2SI+8SJ7s8/X7SyC3LFFXH8sWOr+/VTSk1N7vff7/7Vr37+Grj++riKqybFDPfvtBDuU1vZdrXhDhwJNAANW2yxRTl+DkX1H//hvtlmpdn3c8/Fb+Pcc0uz/3J56624kbXHHpXTtLR8edwUXXvtuNH285/HlVgxLVniPnmy+9e+Fr9HcB861P2cc9znzy/usZqbNi2ON2aMgr0QTU1xlTx8ePzcttnG/cYb3Rsbk66sMGqWKYFcD4TXXy/+vo88MsKnGtqr2zJ5cvycZsxIto6mJvdf/9p9iy2innHj3BctKv1xc0Gf632TC/qf/7z4QX/llbH/b39bwb6mmprc777bfdiw+Blut537rbdWfsgXM9y7AC8DA/JuqO7YyrapDveHH46f2MyZxd3vW29FsB9xRHH3m5QVK9y339594MDkAueZZ9x33z1+Xzvt5P7II8nUsXhxdMfLtffm6jn3XPcXXujYvq+6Kva3337un35anHpr0apV7r/9bTSp5bq/3nFHPF+JihbusS/2AV7I9pr5afa5icCY7Nf/B1gCLAfeBua0tc9qDPf33itN08n558d+K7Gdur1mzoz/02WXlfe477zjfvzxcWN3o43izLZSzsQWLYqrv112+Tzohw1zP+889xdfXLN9XX11/Pt991WwF8uqVe533uk+aJD/60b5b35TeSFf1HAvxUc1hru7+9Zbux90UPH2t2KFe58+7qNGFW+flWL06Oga+eabpT9WY2ME3sYbxw3d445zf/vt0h+3vRYtije+nXf+POgzmXijX7Bg9f92+vTYfp99FOyl0Njofttt0UyTewPu6M33YlK4l8hBB7lvtVXx9nfbbfFbuOee4u2zUsydG2fQxxxT2uP89a8RjOD+jW+4z55d2uMV26uvul966ee9OMC9rs79ggvcX3rpi9tec018f++9i39TWL6osdH9ppvihivEDdh7700+5BXuJXLuufFTe++94uxv553jxVNpl37FcvzxcSb9zDPF3/eSJe7f/a7/a8qAO+5I/g+vo155xf2SS74Y9MOHR9Bfdlk8Hj1awV5OK1dGl8kBA+LnP2KE+x/+kNxrTeFeIrm25Icf7vi+/va32NeUKR3fV6XKdY0cNap4fwyffBLt1D16uK+1lvvPfub+0UfF2XclWbjQ/eKLI0xyQb/XXgr2pKxY4X7ttTEICuLeyYMPlj/kFe4l8sYbXrQbhePGua+3nvsHH3R8X5UsNyT+7rs7tp9c17Wtt/Z/DbF/+eXi1FjpFi50v/lmBXsl+Oyz6KnUt2+8Dr/+9ZjmoFwKDXdN+buGNt0UNtus43O7L1kCd90FRxwRc1Kn2dFHx3KBp57a/tVynn8+FgofOxbWWiumGP7tb2MB81rQv38saL322klXIt26xWRzCxbEhGQvvQS77w7f+hb89a9JV/c5hXs7FGPB7F/8IhYUmDChODVVstyskbk/hjXxwQex2MSQIbF25qRJMePfHnuUplaRQq21Vqy5sGBBvC7nzoXddoM994S//S3p6hTu7VJXB/PmwSeftO/ff/wxXH11nIXWypnn6NFx5j1xYkyz25amJrj+eth223hjOOywmJL3xBNrawpbqXzrrBOvy5dfjiUNZ8+GXXeN1/sTTyRXl8K9HTIZWLUqVn5pj1tugXfeqa4524vh0kvho4/gjDNWv93f/w677AI//CFstVX8gVxzTcytLlKpuneHU06BhQvhwgvhySfhq1+F/faDp54qfz0K93boyILZ7nEJN2wYfP3rxa2r0u2wQ1zGTp8Ozzzz5e+/8Qb84Aex0PHixXDjjfDoo1BfX/5aRdqrR49YkHzhQjjvPPh//y9ew/vvH02K5aJwb4f+/WGDDdoX7n/6U7TNnXRSrBpTa848M352p5zy+YLaK1bEWf2228ZVzU9+AvPnwyGHQCe9QqVK9ewJp58eK3tNnAgPPxwnhgceWJ71efWn0w5m7b+pOnlyNC+MG1f8uqrBRhvFUnUPPQT33AP33x+Lj//4x3Ezas4cuOCC9Pcgktqx3nrws59FyJ9xRvT0amgo/XEV7u2UyUTTQmNj4f/mhRfgvvvgmGPiTnutOuqoaKIZNy5uOjU1wb33xsfAgUlXJ1IaG2wQJzavvBLdWktN4d5OmQx8+mn0vy7UlCnRR/boo0tXVzXo2jUWn+7dO248Pfss7Ltv0lWJlMdGG0GXLqU/ThkOkU75N1UHD257+/feg1/9CsaPh698paSlVYXdd4dXX026CpH00pl7O223XfRvffrpwra/7jpYvjz6w4qIlJrCvZ26dIkbgYXcVG1shKlT44Zh7oxfRKSUFO4dkMlEv9Vcl77W3H13NEHU2qAlEUmOwr0D6urg/fdjsMLqTJ4cfePHjClLWSIiCveOyDWxrK7d/amnYqa444+Hzp3LU5eIiMK9AwYPjsBeXbv75Mmw7rpw+OHlq0tEROHeAWuvDYMGtR7ub7wBt98eMxquv35ZSxORGldQuJvZaDObb2YLzOy0Fr6/lpndkf3+382sf7ELrVR1da2H+5VXwsqVcMIJ5a1JRKTNcDezzsA0YG9gEDDezAY12+xw4F133wa4HLiw2IVWqkwmztBff/2Lz3/6aYT7fvtpSL2IlF8hZ+4jgAXu/rK7rwBuB8Y222YscEP267uAUWa1Medha9P/3n57LEqhQUsikoRCwr0PsDjv8ZLscy1u4+6NwPvAxs13ZGZHmlmDmTUsK2Q5niowbFh8zg/33JztO+4Io0YlU5eI1LZCwr2lM/Dmw3YK2QZ3n+7u9e5e37t370Lqq3jrrQfbbPPF7pCPPBJLbdXqnO0ikrxCwn0J0C/vcV9gaWvbmFkXYH3gnWIUWA2az+0+eTJsvHF5pvUUEWlJIeH+JDDQzAaYWTdgHDCj2TYzgEOzXx8E/Nm9rUH56ZHJxCjV996LRXLvvjvmLF9nnaQrE5Fa1eaUv+7eaGYTgAeAzsAv3X2OmU0EGtx9BnAdcJOZLSDO2GtqnaG6uvg8a1YEe+fOcOyxydYkIrWtoPnc3X0mMLPZc2fkff0p8J3illY9cj1mHnkkpvY9+GDo0/yWs4hIGWkv8rHtAAAE8ElEQVSxjiLYZBPYfHO4+GLN2S4ilUHTDxRJJhPBvssuMGJE0tWISK1TuBdJrt1dc7aLSCVQs0yRHHpozCNzwAFJVyIionAvmq23hvPPT7oKEZGgZhkRkRRSuIuIpJDCXUQkhRTuIiIppHAXEUkhhbuISAop3EVEUkjhLiKSQpbUtOtmtgx4tZ3/vBfwVhHLKRbVtWZU15qr1NpU15rpSF1bunubS9klFu4dYWYN7l6fdB3Nqa41o7rWXKXWprrWTDnqUrOMiEgKKdxFRFKoWsN9etIFtEJ1rRnVteYqtTbVtWZKXldVtrmLiMjqVeuZu4iIrEZVhbuZ9TOzv5jZPDObY2YVsVqpma1tZk+Y2exsXWcnXVM+M+tsZv8ws3uTriXHzF4xs2fNbJaZNSRdT46ZbWBmd5nZ89nX2S4VUNN22Z9T7uMDM6uINb/M7OTsa/45M7vNzNZOuiYAMzsxW9OcJH9WZvZLM3vTzJ7Le24jM3vQzF7Mft6wFMeuqnAHGoFT3X0HYGfgODMblHBNAJ8B33L3nYBhwGgz2znhmvKdCMxLuogW7O7uwyqsq9pk4H533x7YiQr4ubn7/OzPaRgwHPgY+F3CZWFmfYATgHp3Hwx0BsYlWxWY2WDgR8AI4ne4n5kNTKicXwGjmz13GvCQuw8EHso+LrqqCnd3f93dn85+/SHxh9cn2arAw0fZh12zHxVxM8PM+gL7AtcmXUulM7P1gN2A6wDcfYW7v5dsVV8yCnjJ3ds7ALDYugDrmFkXoDuwNOF6AHYAHnf3j929EXgESGQBTHf/X+CdZk+PBW7Ifn0DsH8pjl1V4Z7PzPoDGeDvyVYSsk0fs4A3gQfdvSLqAiYB/wU0JV1IMw780cyeMrMjky4maytgGXB9thnrWjPrkXRRzYwDbku6CAB3fw24BFgEvA687+5/TLYqAJ4DdjOzjc2sO7AP0C/hmvJt6u6vQ5ywApuU4iBVGe5mti7wG+Akd/8g6XoA3H1V9rK5LzAie2mYKDPbD3jT3Z9KupYWjHT3OmBvonltt6QLIs5C64Ar3T0DLKdEl8ztYWbdgDHAr5OuBSDbVjwWGABsDvQws+8nWxW4+zzgQuBB4H5gNtGkW1OqLtzNrCsR7Le4+2+Trqe57GX8w3y5nS0JI4ExZvYKcDvwLTO7OdmSgrsvzX5+k2g/HpFsRQAsAZbkXXXdRYR9pdgbeNrd/5l0IVl7AAvdfZm7rwR+C+yacE0AuPt17l7n7rsRzSIvJl1Tnn+a2WYA2c9vluIgVRXuZmZEe+g8d78s6XpyzKy3mW2Q/Xod4kX/fLJVgbuf7u593b0/cTn/Z3dP/MzKzHqYWc/c18C/EZfSiXL3N4DFZrZd9qlRwNwES2puPBXSJJO1CNjZzLpn/zZHUQE3oAHMbJPs5y2A/0tl/dxmAIdmvz4UuLsUB+lSip2W0EjgEODZbPs2wH+7+8wEawLYDLjBzDoTb5h3unvFdDusQJsCv4s8oAtwq7vfn2xJ/3I8cEu2CeRl4AcJ1wNAtu14T+CopGvJcfe/m9ldwNNEs8c/qJwRob8xs42BlcBx7v5uEkWY2W3AN4FeZrYEOBO4ALjTzA4n3iC/U5Jja4SqiEj6VFWzjIiIFEbhLiKSQgp3EZEUUriLiKSQwl1EJIUU7iIiKaRwFxFJIYW7iEgK/X/Ufqnf+csVJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x33009c42e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 绘制不同聚类数目的模型性能， 找到最佳模型\n",
    "plt.plot(Ks, np.array(CH_scores), 'b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "K为1会报错，K为2时分数最好，不太明白。。0.5以上的取值只有K为2这一种。。是简单分成室内活动和室外活动了么。。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最近加班严重，本次作业没有时间仔细看了。。掉队了。。。难受"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
